<!DOCTYPE html>
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<meta name="generator" content="hevea 2.09">

<META name="Author" content="Julien Mairal">
<link rel="stylesheet" href="doc_spams.css"><link rel="stylesheet" type="text/css" href="doc_spams.css">
<title>Proximal Toolbox</title>
</head>
<body>
<a href="doc_spams005.html"><img src="previous_motif.gif" alt="Previous"></a>
<a href="index.html"><img src="contents_motif.gif" alt="Up"></a>
<a href="doc_spams007.html"><img src="next_motif.gif" alt="Next"></a>
<hr>
<h2 class="section" id="sec23">5  Proximal Toolbox</h2>
<ul>
<li><a href="doc_spams006.html#sec24">Regularization Functions</a>
</li><li><a href="doc_spams006.html#sec25">Function mexProximalFlat</a>
</li><li><a href="doc_spams006.html#sec26">Function mexProximalTree</a>
</li><li><a href="doc_spams006.html#sec27">Function mexProximalGraph</a>
</li><li><a href="doc_spams006.html#sec28">Function mexProximalPathCoding</a>
</li><li><a href="doc_spams006.html#sec29">Function mexEvalPathCoding</a>
</li><li><a href="doc_spams006.html#sec30">Problems Addressed</a>
</li><li><a href="doc_spams006.html#sec37">Function mexFistaFlat</a>
</li><li><a href="doc_spams006.html#sec38">Function mexFistaTree</a>
</li><li><a href="doc_spams006.html#sec39">Function mexFistaGraph</a>
</li><li><a href="doc_spams006.html#sec40">Function mexFistaPathCoding</a>
</li><li><a href="doc_spams006.html#sec41">Function solverPoisson</a>
</li><li><a href="doc_spams006.html#sec42">Function mexIncrementalProx</a>
</li><li><a href="doc_spams006.html#sec43">Function mexStochasticProx</a>
</li></ul>
<p>
The previous toolbox we have presented is well
adapted for solving a large number of small and medium-scale sparse
decomposition problems with the square loss, which is typical from the
classical dictionary learning framework. We now present
a new software package that is adapted for solving a wide range of
possibly large-scale learning problems, with several combinations of losses and
regularization terms. The method implements the proximal methods
of [<a href="doc_spams010.html#beck">1</a>], and includes the proximal solvers for the tree-structured
regularization of [<a href="doc_spams010.html#jenatton3">15</a>], and the solver of [<a href="doc_spams010.html#mairal10">22</a>] for
general structured sparse regularization. 
The solver for structured sparse regularization norms includes a C++ max-flow
implementation of the push-relabel algorithm of [<a href="doc_spams010.html#goldberg">13</a>], with
heuristics proposed by [<a href="doc_spams010.html#cherkassky">5</a>].</p><p>This implementation also provides robust stopping criteria based on
<em>duality gaps</em>, which are presented in Appendix <a href="doc_spams009.html#appendix">A</a>. It can handle intercepts (unregularized variables). The general formulation that our software
can solve take the form
</p><table class="display dcenter"><tr class="c016"><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">w</span> ∈ ℝ<sup><span class="c007">p</span></sup></td></tr>
</table></td><td class="dcell"> [<span class="c007">g</span>(<span class="c009">w</span>) </td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c005">▵</span></td></tr>
<tr><td class="dcell c012">=</td></tr>
<tr><td class="dcell c012">&nbsp;</td></tr>
</table></td><td class="dcell"> <span class="c007">f</span>(<span class="c009">w</span>) + λψ(<span class="c009">w</span>)],
</td></tr>
</table><p>
where <span class="c007">f</span> is a smooth loss function and ψ is a regularization function.
When one optimizes a matrix <span class="c009">W</span> in ℝ<sup><span class="c007">p</span> × <span class="c007">r</span></sup> instead of
a vector <span class="c009">w</span> in ℝ<sup><span class="c007">p</span></sup>, we will write 
</p><table class="display dcenter"><tr class="c016"><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">W</span> ∈ ℝ<sup><span class="c007">p</span> × <span class="c007">r</span></sup></td></tr>
</table></td><td class="dcell"> [<span class="c007">g</span>(<span class="c009">W</span>) </td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c005">▵</span></td></tr>
<tr><td class="dcell c012">=</td></tr>
<tr><td class="dcell c012">&nbsp;</td></tr>
</table></td><td class="dcell"> <span class="c007">f</span>(<span class="c009">W</span>) + λψ(<span class="c009">W</span>)].
</td></tr>
</table><p>
Note that the software can possibly handle nonnegativity constraints.</p><p>We start by presenting the type of regularization implemented in the software
</p>
<h3 class="subsection" id="sec24">5.1  Regularization Functions</h3>
<p>
Our software can handle the following regularization functions ψ for vectors <span class="c009">w</span> in ℝ<sup><span class="c007">p</span></sup>:
</p><ul class="itemize"><li class="li-itemize">
<span class="c010">The Tikhonov regularization</span>: ψ(<span class="c009">w</span>) =<sup><span class="c005">▵</span></sup> 1/2||<span class="c009">w</span>||<sub>2</sub><sup>2</sup>.
</li><li class="li-itemize"><span class="c010">The </span>ℓ<sub>1</sub><span class="c010">-norm</span>: ψ(<span class="c009">w</span>) =<sup><span class="c005">▵</span></sup> ||<span class="c009">w</span>||<sub>1</sub>.
</li><li class="li-itemize"><span class="c010">The Elastic-Net</span>: ψ(<span class="c009">w</span>) =<sup><span class="c005">▵</span></sup> ||<span class="c009">w</span>||<sub>1</sub>+γ||<span class="c009">w</span>||<sub>2</sub><sup>2</sup>.
</li><li class="li-itemize"><span class="c010">The Fused-Lasso</span>: ψ(<span class="c009">w</span>) =<sup><span class="c005">▵</span></sup> ||<span class="c009">w</span>||<sub>1</sub>+γ||<span class="c009">w</span>||<sub>2</sub><sup>2</sup>+γ<sub>2</sub>∑<sub><span class="c007">i</span>=1</sub><sup><span class="c007">p</span>−1</sup>|<span class="c009">w</span><sub><span class="c007">i</span>+1</sub>−<span class="c009">w</span><sub><span class="c007">i</span></sub>|.
</li><li class="li-itemize"><span class="c010">The group Lasso</span>: ψ(<span class="c009">w</span>) =<sup><span class="c005">▵</span></sup> ∑<sub><span class="c007">g</span> ∈ <span class="c008">G</span></sub> η<sub><span class="c007">g</span></sub> ||<span class="c009">w</span><sub><span class="c007">g</span></sub>||<sub>2</sub>, where <span class="c008">G</span> are groups of variables.
</li><li class="li-itemize"><span class="c010">The group Lasso with </span>ℓ<sub>∞</sub><span class="c010">-norm</span>: ψ(<span class="c009">w</span>) =<sup><span class="c005">▵</span></sup> ∑<sub><span class="c007">g</span> ∈ <span class="c008">G</span></sub> η<sub><span class="c007">g</span></sub> ||<span class="c009">w</span><sub><span class="c007">g</span></sub>||<sub>∞</sub>, where <span class="c008">G</span> are groups of variables.
</li><li class="li-itemize"><span class="c010">The sparse group Lasso</span>: same as above but with an additional ℓ<sub>1</sub> term.
</li><li class="li-itemize"><span class="c010">The tree-structured sum of </span>ℓ<sub>2</sub><span class="c010">-norms</span>: ψ(<span class="c009">w</span>) =<sup><span class="c005">▵</span></sup> ∑<sub><span class="c007">g</span> ∈ <span class="c008">G</span></sub> η<sub><span class="c007">g</span></sub> ||<span class="c009">w</span><sub><span class="c007">g</span></sub>||<sub>2</sub>, where <span class="c008">G</span> is a tree-structured set of groups [<a href="doc_spams010.html#jenatton3">15</a>], and the η<sub><span class="c007">g</span></sub> are positive weights.
</li><li class="li-itemize"><span class="c010">The tree-structured sum of </span>ℓ<sub>∞</sub><span class="c010">-norms</span>: ψ(<span class="c009">w</span>) =<sup><span class="c005">▵</span></sup> ∑<sub><span class="c007">g</span> ∈ <span class="c008">G</span></sub> η<sub><span class="c007">g</span></sub> ||<span class="c009">w</span><sub><span class="c007">g</span></sub>||<sub>∞</sub>. See [<a href="doc_spams010.html#jenatton3">15</a>]
</li><li class="li-itemize"><span class="c010">General sum of </span>ℓ<sub>∞</sub><span class="c010">-norms</span>: ψ(<span class="c009">w</span>) =<sup><span class="c005">▵</span></sup> ∑<sub><span class="c007">g</span> ∈ <span class="c008">G</span></sub> η<sub><span class="c007">g</span></sub> ||<span class="c009">w</span><sub><span class="c007">g</span></sub>||<sub>∞</sub>, where no assumption are made on the groups <span class="c008">G</span>.
</li><li class="li-itemize"><span class="c010">The path-coding penalties</span> of [<a href="doc_spams010.html#mairal14">24</a>].
</li><li class="li-itemize"><span class="c010">the </span>ℓ<sub>1</sub><span class="c010">-constraint</span>.
</li><li class="li-itemize"><span class="c010">a few undocumented ones</span>.
</li></ul><p>
Our software also handles regularization functions ψ on matrices <span class="c009">W</span> in ℝ<sup><span class="c007">p</span> × <span class="c007">r</span></sup> (note that <span class="c009">W</span> can be transposed in these formulations). In particular,
</p><ul class="itemize"><li class="li-itemize">
<span class="c010">The </span>ℓ<sub>1</sub>/ℓ<sub>2</sub><span class="c010">-norm</span>: ψ(<span class="c009">W</span>) =<sup><span class="c005">▵</span></sup> ∑<sub><span class="c007">i</span>=1</sub><sup><span class="c007">p</span></sup> ||<span class="c009">W</span><sub><span class="c007">i</span></sub>||<sub>2</sub>, where <span class="c009">W</span><sub><span class="c007">i</span></sub> denotes the <span class="c007">i</span>-th row of <span class="c009">W</span>.
</li><li class="li-itemize"><span class="c010">The </span>ℓ<sub>1</sub>/ℓ<sub>∞</sub><span class="c010">-norm</span>: ψ(<span class="c009">W</span>) =<sup><span class="c005">▵</span></sup> ∑<sub><span class="c007">i</span>=1</sub><sup><span class="c007">p</span></sup> ||<span class="c009">W</span><sub><span class="c007">i</span></sub>||<sub>∞</sub>,
</li><li class="li-itemize"><span class="c010">The </span>ℓ<sub>1</sub>/ℓ<sub>2</sub><span class="c010">+</span>ℓ<sub>1</sub><span class="c010">-norm</span>: ψ(<span class="c009">W</span>) =<sup><span class="c005">▵</span></sup> ∑<sub><span class="c007">i</span>=1</sub><sup><span class="c007">p</span></sup> ||<span class="c009">W</span><sub><span class="c007">i</span></sub>||<sub>2</sub> + λ<sub>2</sub> ∑<sub><span class="c007">i</span>,<span class="c007">j</span></sub>|<span class="c009">W</span><sub><span class="c007">ij</span></sub>|.
</li><li class="li-itemize"><span class="c010">The </span>ℓ<sub>1</sub>/ℓ<sub>∞</sub><span class="c010">+</span>ℓ<sub>1</sub><span class="c010">-norm</span>: ψ(<span class="c009">W</span>) =<sup><span class="c005">▵</span></sup> ∑<sub><span class="c007">i</span>=1</sub><sup><span class="c007">p</span></sup> ||<span class="c009">W</span><sub><span class="c007">i</span></sub>||<sub>∞</sub>+λ<sub>2</sub> ∑<sub><span class="c007">i</span>,<span class="c007">j</span></sub>|<span class="c009">W</span><sub><span class="c007">ij</span></sub>|,
</li><li class="li-itemize"><span class="c010">The </span>ℓ<sub>1</sub>/ℓ<sub>∞</sub><span class="c010">-norm on rows and columns</span>: ψ(<span class="c009">W</span>) =<sup><span class="c005">▵</span></sup> ∑<sub><span class="c007">i</span>=1</sub><sup><span class="c007">p</span></sup> ||<span class="c009">W</span><sub><span class="c007">i</span></sub>||<sub>∞</sub>+λ<sub>2</sub> ∑<sub><span class="c007">j</span>=1</sub><sup><span class="c007">r</span></sup>||<span class="c009">W</span><sup><span class="c007">j</span></sup>||<sub>∞</sub>, where <span class="c009">W</span><sup><span class="c007">j</span></sup> denotes the <span class="c007">j</span>-th column of <span class="c009">W</span>.
</li><li class="li-itemize"><span class="c010">The multi-task tree-structured sum of </span>ℓ<sub>∞</sub><span class="c010">-norms</span>: 
<table class="display dcenter"><tr class="c016"><td class="dcell">
ψ(<span class="c009">W</span>)</td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c005">▵</span></td></tr>
<tr><td class="dcell c012">=</td></tr>
<tr><td class="dcell c012">&nbsp;</td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c007">r</span></td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">i</span>=1</td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">g</span> ∈ <span class="c008">G</span></td></tr>
</table></td><td class="dcell"> η<sub><span class="c007">g</span></sub>||<span class="c009">w</span><sub><span class="c007">g</span></sub><sup><span class="c007">i</span></sup>||<sub>∞</sub>+ γ </td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">g</span> ∈ <span class="c008">G</span></td></tr>
</table></td><td class="dcell"> η<sub><span class="c007">g</span></sub> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">max</td></tr>
<tr><td class="dcell c012"><span class="c007">j</span> ∈ <span class="c007">g</span></td></tr>
</table></td><td class="dcell">||<span class="c009">W</span><sub><span class="c007">j</span></sub>||<sub>∞</sub>,
<a id="software:eq:struct"></a>
    (25)</td></tr>
</table>
where the first double sums is in fact a sum of independent structured norms on the columns <span class="c009">w</span><sup><span class="c007">i</span></sup> of <span class="c009">W</span>, and the right term is a tree-structured regularization norm applied to the ℓ<sub>∞</sub>-norm of the rows of <span class="c009">W</span>, thereby inducing the tree-structured regularization at the row level. <span class="c008">G</span> is here a tree-structured set of groups.
</li><li class="li-itemize"><span class="c010">The multi-task general sum of </span>ℓ<sub>∞</sub><span class="c010">-norms</span> is the same as Eq. (<a href="#software%3Aeq%3Astruct">25</a>) except that the groups <span class="c008">G</span> are general overlapping groups.
</li><li class="li-itemize"><span class="c010">The trace norm</span>: ψ(<span class="c009">W</span>) =<sup><span class="c005">▵</span></sup> ||<span class="c009">W</span>||<sub>*</sub>.
</li></ul><p>
Non-convex regularizations are also implemented with the ISTA algorithm (no duality gaps are of course provided in these cases):
</p><ul class="itemize"><li class="li-itemize">
<span class="c010">The </span>ℓ<sub>0</sub><span class="c010">-pseudo-norm</span>: ψ(<span class="c009">w</span>) =<sup><span class="c005">▵</span></sup> ||<span class="c009">w</span>||<sub>0</sub>.
</li><li class="li-itemize"><span class="c010">The rank</span>: ψ(<span class="c009">W</span>) =<sup><span class="c005">▵</span></sup> randk(<span class="c009">W</span>).
</li><li class="li-itemize"><span class="c010">The tree-structured </span>ℓ<sub>0</sub><span class="c010">-pseudo-norm</span>: ψ(<span class="c009">w</span>) =<sup><span class="c005">▵</span></sup> ∑<sub><span class="c007">g</span> ∈ <span class="c008">G</span></sub> δ<sub><span class="c009">w</span><sub><span class="c007">g</span></sub> ≠ 0</sub>.
</li></ul><p>All of these regularization terms for vectors or matrices can be coupled with
nonnegativity constraints. It is also possible to add an intercept, which one
wishes not to regularize, and we will include this possibility in the next
sections. There are also a few hidden undocumented options which are available in the source code.</p><p>We now present 3 functions for computing proximal operators associated to the previous regularization functions.
</p>
<h3 class="subsection" id="sec25">5.2  Function mexProximalFlat</h3>
<p>
This function computes the proximal operators associated to many regularization functions, for input signals <span class="c009">U</span>=[<span class="c009">u</span><sup>1</sup>,…,<span class="c009">u</span><sup><span class="c007">n</span></sup>] in ℝ<sup><span class="c007">p</span> × <span class="c007">n</span></sup>, it finds a matrix <span class="c009">V</span>=[<span class="c009">v</span><sup>1</sup>,…,<span class="c009">v</span><sup><span class="c007">n</span></sup>] in ℝ<sup><span class="c007">p</span> × <span class="c007">n</span></sup> such that:</p><p>•  If one chooses a regularization function on vectors, for every column <span class="c009">u</span> of <span class="c009">U</span>, it computes one column <span class="c009">v</span> of <span class="c009">V</span> solving
</p><table class="display dcenter"><tr class="c016"><td class="dcell">
</td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">v</span> ∈ ℝ<sup><span class="c007">p</span></sup></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012">2</td></tr>
</table></td><td class="dcell">||<span class="c009">u</span>−<span class="c009">v</span>||<sub>2</sub><sup>2</sup> + λ ||<span class="c009">v</span>||<sub>0</sub>,
    (26)</td></tr>
</table><p>
or
</p><table class="display dcenter"><tr class="c016"><td class="dcell">
</td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">v</span> ∈ ℝ<sup><span class="c007">p</span></sup></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012">2</td></tr>
</table></td><td class="dcell">||<span class="c009">u</span>−<span class="c009">v</span>||<sub>2</sub><sup>2</sup> + λ ||<span class="c009">v</span>||<sub>1</sub>,
    (27)</td></tr>
</table><p>
or
</p><table class="display dcenter"><tr class="c016"><td class="dcell">
</td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">v</span> ∈ ℝ<sup><span class="c007">p</span></sup></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012">2</td></tr>
</table></td><td class="dcell">||<span class="c009">u</span>−<span class="c009">v</span>||<sub>2</sub><sup>2</sup> + </td><td class="dcell"><table class="display"><tr><td class="dcell c012">λ</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012">2</td></tr>
</table></td><td class="dcell"> ||<span class="c009">v</span>||<sub>2</sub><sup>2</sup>,
    (28)</td></tr>
</table><p>
or
</p><table class="display dcenter"><tr class="c016"><td class="dcell">
</td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">v</span> ∈ ℝ<sup><span class="c007">p</span></sup></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012">2</td></tr>
</table></td><td class="dcell">||<span class="c009">u</span>−<span class="c009">v</span>||<sub>2</sub><sup>2</sup> + λ ||<span class="c009">v</span>||<sub>1</sub> + λ<sub>2</sub>||<span class="c009">v</span>||<sub>2</sub><sup>2</sup>,
    (29)</td></tr>
</table><p>
or
</p><table class="display dcenter"><tr class="c016"><td class="dcell">
</td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">v</span> ∈ ℝ<sup><span class="c007">p</span></sup></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012">2</td></tr>
</table></td><td class="dcell">||<span class="c009">u</span>−<span class="c009">v</span>||<sub>2</sub><sup>2</sup> + λ</td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c007">p</span>−1</td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">j</span>=1</td></tr>
</table></td><td class="dcell">|<span class="c009">v</span><sub><span class="c007">j</span>+1</sub><sup><span class="c007">i</span></sup>−<span class="c009">v</span><sub><span class="c007">j</span></sub><sup><span class="c007">i</span></sup>|+λ<sub>2</sub> ||<span class="c009">v</span>||<sub>1</sub> + λ<sub>3</sub>||<span class="c009">v</span>||<sub>2</sub><sup>2</sup>,
    (30)</td></tr>
</table><p>
or
</p><table class="display dcenter"><tr class="c016"><td class="dcell">
</td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">v</span> ∈ ℝ<sup><span class="c007">p</span></sup></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012">2</td></tr>
</table></td><td class="dcell">||<span class="c009">u</span>−<span class="c009">v</span>||<sub>2</sub><sup>2</sup> + λ </td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">g</span> ∈ <span class="c008">T</span></td></tr>
</table></td><td class="dcell"> δ<sup><span class="c007">g</span></sup>(<span class="c009">v</span>),
    (31)</td></tr>
</table><p>
where <span class="c008">T</span> is a tree-structured set of groups (see [<a href="doc_spams010.html#jenatton4">16</a>]), and δ<sup><span class="c007">g</span></sup>(<span class="c009">v</span>) = 0 if <span class="c009">v</span><sub><span class="c007">g</span></sub>=0 and 1 otherwise.
It can also solve
</p><table class="display dcenter"><tr class="c016"><td class="dcell">
</td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">v</span> ∈ ℝ<sup><span class="c007">p</span></sup></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012">2</td></tr>
</table></td><td class="dcell">||<span class="c009">u</span>−<span class="c009">v</span>||<sub>2</sub><sup>2</sup> + λ </td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">g</span> ∈ <span class="c008">T</span></td></tr>
</table></td><td class="dcell"> η<sup><span class="c007">g</span></sup> ||<span class="c009">v</span><sub><span class="c007">g</span></sub>||<sub>2</sub>,
    (32)</td></tr>
</table><p>
or
</p><table class="display dcenter"><tr class="c016"><td class="dcell">
</td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">v</span> ∈ ℝ<sup><span class="c007">p</span></sup></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012">2</td></tr>
</table></td><td class="dcell">||<span class="c009">u</span>−<span class="c009">v</span>||<sub>2</sub><sup>2</sup> + λ </td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">g</span> ∈ <span class="c008">T</span></td></tr>
</table></td><td class="dcell"> η<sup><span class="c007">g</span></sup> ||<span class="c009">v</span><sub><span class="c007">g</span></sub>||<sub>∞</sub>,
    (33)</td></tr>
</table><p>
or
</p><table class="display dcenter"><tr class="c016"><td class="dcell">
</td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">v</span> ∈ ℝ<sup><span class="c007">p</span></sup></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012">2</td></tr>
</table></td><td class="dcell">||<span class="c009">u</span>−<span class="c009">v</span>||<sub>2</sub><sup>2</sup> + λ </td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">g</span> ∈ <span class="c008">G</span></td></tr>
</table></td><td class="dcell"> η<sup><span class="c007">g</span></sup> ||<span class="c009">v</span><sub><span class="c007">g</span></sub>||<sub>∞</sub>,
    (34)</td></tr>
</table><p>
where <span class="c008">G</span> is any kind of set of groups.</p><p>This function can also solve the following proximal operators on matrices
</p><table class="display dcenter"><tr class="c016"><td class="dcell">
</td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">V</span> ∈ ℝ<sup><span class="c007">p</span> × <span class="c007">n</span></sup></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012">2</td></tr>
</table></td><td class="dcell">||<span class="c009">U</span>−<span class="c009">V</span>||<sub><span class="c007">F</span></sub><sup>2</sup> + λ </td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c007">p</span></td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">i</span>=1</td></tr>
</table></td><td class="dcell"> ||<span class="c009">V</span><sub><span class="c007">i</span></sub>||<sub>2</sub>, 
    (35)</td></tr>
</table><p>
where <span class="c009">V</span><sub><span class="c007">i</span></sub> is the <span class="c007">i</span>-th row of <span class="c009">V</span>, or
</p><table class="display dcenter"><tr class="c016"><td class="dcell">
</td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">V</span> ∈ ℝ<sup><span class="c007">p</span> × <span class="c007">n</span></sup></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012">2</td></tr>
</table></td><td class="dcell">||<span class="c009">U</span>−<span class="c009">V</span>||<sub><span class="c007">F</span></sub><sup>2</sup> + λ </td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c007">p</span></td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">i</span>=1</td></tr>
</table></td><td class="dcell"> ||<span class="c009">V</span><sub><span class="c007">i</span></sub>||<sub>∞</sub>, 
    (36)</td></tr>
</table><p>
or
</p><table class="display dcenter"><tr class="c016"><td class="dcell">
</td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">V</span> ∈ ℝ<sup><span class="c007">p</span> × <span class="c007">n</span></sup></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012">2</td></tr>
</table></td><td class="dcell">||<span class="c009">U</span>−<span class="c009">V</span>||<sub><span class="c007">F</span></sub><sup>2</sup> + λ </td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c007">p</span></td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">i</span>=1</td></tr>
</table></td><td class="dcell"> ||<span class="c009">V</span><sub><span class="c007">i</span></sub>||<sub>2</sub> +λ<sub>2</sub> </td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c007">p</span></td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">i</span>=1</td></tr>
</table></td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c007">n</span></td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">j</span>=1</td></tr>
</table></td><td class="dcell"> |<span class="c009">V</span><sub><span class="c007">ij</span></sub>|, 
    (37)</td></tr>
</table><p>
or
</p><table class="display dcenter"><tr class="c016"><td class="dcell">
</td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">V</span> ∈ ℝ<sup><span class="c007">p</span> × <span class="c007">n</span></sup></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012">2</td></tr>
</table></td><td class="dcell">||<span class="c009">U</span>−<span class="c009">V</span>||<sub><span class="c007">F</span></sub><sup>2</sup> + λ </td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c007">p</span></td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">i</span>=1</td></tr>
</table></td><td class="dcell"> ||<span class="c009">V</span><sub><span class="c007">i</span></sub>||<sub>∞</sub>+λ<sub>2</sub> </td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c007">p</span></td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">i</span>=1</td></tr>
</table></td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c007">n</span></td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">j</span>=1</td></tr>
</table></td><td class="dcell"> |<span class="c009">V</span><sub><span class="c007">ij</span></sub>|, 
    (38)</td></tr>
</table><p>
or
</p><table class="display dcenter"><tr class="c016"><td class="dcell">
</td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">V</span> ∈ ℝ<sup><span class="c007">p</span> × <span class="c007">n</span></sup></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012">2</td></tr>
</table></td><td class="dcell">||<span class="c009">U</span>−<span class="c009">V</span>||<sub><span class="c007">F</span></sub><sup>2</sup> + λ </td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c007">p</span></td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">i</span>=1</td></tr>
</table></td><td class="dcell"> ||<span class="c009">V</span><sub><span class="c007">i</span></sub>||<sub>∞</sub>+λ<sub>2</sub> </td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c007">n</span></td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">j</span>=1</td></tr>
</table></td><td class="dcell"> ||<span class="c009">V</span><sup><span class="c007">j</span></sup>||<sub>∞</sub>.
    (39)</td></tr>
</table><p>
where <span class="c009">V</span><sup><span class="c007">j</span></sup> is the <span class="c007">j</span>-th column of <span class="c009">V</span>.</p><p>See details below:
</p><table class="lstframe c011"><tr><td class="mouselstlisting"><span class="c001">% <br>
% Usage:  [V [val_regularizer]]=mexProximalFlat(U,param);<br>
%<br>
% Name: mexProximalFlat<br>
%<br>
% Description: mexProximalFlat computes proximal operators. Depending<br>
%         on the value of param.regul, it computes <br>
%<br>
%         Given an input matrix U=[u^1,\ldots,u^n], it computes a matrix <br>
%         V=[v^1,\ldots,v^n] such that<br>
%         if one chooses a regularization functions on vectors, it computes<br>
%         for each column u of U, a column v of V solving<br>
%         if param.regul='l0'<br>
%             argmin 0.5||u-v||_2^2 + lambda||v||_0<br>
%         if param.regul='l1'<br>
%             argmin 0.5||u-v||_2^2 + lambda||v||_1<br>
%         if param.regul='l2'<br>
%             argmin 0.5||u-v||_2^2 + 0.5lambda||v||_2^2<br>
%         if param.regul='elastic-net'<br>
%             argmin 0.5||u-v||_2^2 + lambda||v||_1 + lambda_2||v||_2^2<br>
%         if param.regul='fused-lasso'<br>
%             argmin 0.5||u-v||_2^2 + lambda FL(v) + ...<br>
%                               ...  lambda_2||v||_1 + lambda_3||v||_2^2<br>
%         if param.regul='linf'<br>
%             argmin 0.5||u-v||_2^2 + lambda||v||_inf<br>
%         if param.regul='l1-constraint'<br>
%             argmin 0.5||u-v||_2^2 s.t. ||v||_1 &lt;= lambda<br>
%         if param.regul='l2-not-squared'<br>
%             argmin 0.5||u-v||_2^2 + lambda||v||_2<br>
%         if param.regul='group-lasso-l2'  <br>
%             argmin 0.5||u-v||_2^2 + lambda sum_g ||v_g||_2 <br>
%             where the groups are either defined by param.groups or by param.size_group,<br>
%         if param.regul='group-lasso-linf'<br>
%             argmin 0.5||u-v||_2^2 + lambda sum_g ||v_g||_inf<br>
%         if param.regul='sparse-group-lasso-l2'  <br>
%             argmin 0.5||u-v||_2^2 + lambda sum_g ||v_g||_2 + lambda_2 ||v||_1<br>
%             where the groups are either defined by param.groups or by param.size_group,<br>
%         if param.regul='sparse-group-lasso-linf'<br>
%             argmin 0.5||u-v||_2^2 + lambda sum_g ||v_g||_inf + lambda_2 ||v||_1<br>
%         if param.regul='trace-norm-vec' <br>
%             argmin 0.5||u-v||_2^2 + lambda ||mat(v)||_* <br>
%            where mat(v) has param.size_group rows<br>
%<br>
%         if one chooses a regularization function on matrices<br>
%         if param.regul='l1l2',  V= <br>
%             argmin 0.5||U-V||_F^2 + lambda||V||_{1/2}<br>
%         if param.regul='l1linf',  V= <br>
%             argmin 0.5||U-V||_F^2 + lambda||V||_{1/inf}<br>
%         if param.regul='l1l2+l1',  V= <br>
%             argmin 0.5||U-V||_F^2 + lambda||V||_{1/2} + lambda_2||V||_{1/1}<br>
%         if param.regul='l1linf+l1',  V= <br>
%             argmin 0.5||U-V||_F^2 + lambda||V||_{1/inf} + lambda_2||V||_{1/1}<br>
%         if param.regul='l1linf+row-column',  V= <br>
%             argmin 0.5||U-V||_F^2 + lambda||V||_{1/inf} + lambda_2||V'||_{1/inf}<br>
%         if param.regul='trace-norm',  V= <br>
%             argmin 0.5||U-V||_F^2 + lambda||V||_*<br>
%         if param.regul='rank',  V= <br>
%             argmin 0.5||U-V||_F^2 + lambda rank(V)<br>
%         if param.regul='none',  V= <br>
%             argmin 0.5||U-V||_F^2 <br>
%         <br>
%         for all these regularizations, it is possible to enforce non-negativity constraints<br>
%         with the option param.pos, and to prevent the last row of U to be regularized, with<br>
%         the option param.intercept<br>
%<br>
% Inputs: U:  double m x n matrix   (input signals)<br>
%               m is the signal size<br>
%         param: struct<br>
%               param.lambda  (regularization parameter)<br>
%               param.regul (choice of regularization, see above)<br>
%               param.lambda2  (optional, regularization parameter)<br>
%               param.lambda3  (optional, regularization parameter)<br>
%               param.verbose (optional, verbosity level, false by default)<br>
%               param.intercept (optional, last row of U is not regularized,<br>
%                 false by default)<br>
%               param.transpose (optional, transpose the matrix in the regularization function)<br>
%               param.size_group (optional, for regularization functions assuming a group<br>
%                 structure). It is a scalar. When param.groups is not specified, it assumes<br>
%                 that the groups are the sets of consecutive elements of size param.size_group<br>
%               param.groups (int32, optional, for regularization functions assuming a group<br>
%                 structure. It is an int32 vector of size m containing the group indices of the<br>
%                 variables (first group is 1).<br>
%               param.pos (optional, adds positivity constraints on the<br>
%                 coefficients, false by default)<br>
%               param.numThreads (optional, number of threads for exploiting<br>
%                 multi-core / multi-cpus. By default, it takes the value -1,<br>
%                 which automatically selects all the available CPUs/cores).<br>
%<br>
% Output: V: double m x n matrix (output coefficients)<br>
%         val_regularizer: double 1 x n vector (value of the regularization<br>
%         term at the optimum).<br>
%<br>
% Author: Julien Mairal, 2010</span></td></tr>
</table><p>
The following piece of code illustrates how to use this function.
</p><table class="lstframe c011"><tr><td class="mouselstlisting">U=<span class="c002">randn</span>(100,1000);<br>
<br>
param.lambda=0.1; <span class="c001">% regularization parameter</span><br>
param.num_threads=-1; <span class="c001">% all cores (-1 by default)</span><br>
param.verbose=true;   <span class="c001">% verbosity, false by default<br>
<br>
% test l0</span><br>
<span class="c002">fprintf</span>(<span class="c003">'\nprox l0\n'</span>);<br>
param.regul=<span class="c003">'l0'</span>;<br>
param.pos=false;       <span class="c001">% false by default</span><br>
param.intercept=false; <span class="c001">% false by default</span><br>
alpha=mexProximalFlat(U,param);<br>
<br>
<span class="c001">% test l1</span><br>
<span class="c002">fprintf</span>(<span class="c003">'\nprox l1, intercept, positivity constraint\n'</span>);<br>
param.regul=<span class="c003">'l1'</span>;<br>
param.pos=true;       <span class="c001">% can be used with all the other regularizations</span><br>
param.intercept=true; <span class="c001">% can be used with all the other regularizations     </span><br>
alpha=mexProximalFlat(U,param);<br>
<br>
<span class="c001">% test l2</span><br>
<span class="c002">fprintf</span>(<span class="c003">'\nprox squared-l2\n'</span>);<br>
param.regul=<span class="c003">'l2'</span>;<br>
param.pos=false;<br>
param.intercept=false;<br>
alpha=mexProximalFlat(U,param);<br>
<br>
<span class="c001">% test elastic-net</span><br>
<span class="c002">fprintf</span>(<span class="c003">'\nprox elastic-net\n'</span>);<br>
param.regul=<span class="c003">'elastic-net'</span>;<br>
param.lambda2=0.1;<br>
alpha=mexProximalFlat(U,param);<br>
<br>
<span class="c001">% test fused-lasso</span><br>
<span class="c002">fprintf</span>(<span class="c003">'\nprox fused lasso\n'</span>);<br>
param.regul=<span class="c003">'fused-lasso'</span>;<br>
param.lambda2=0.1;<br>
param.lambda3=0.1;<br>
alpha=mexProximalFlat(U,param);<br>
<br>
<span class="c001">% test l1l2</span><br>
<span class="c002">fprintf</span>(<span class="c003">'\nprox mixed norm l1/l2\n'</span>);<br>
param.regul=<span class="c003">'l1l2'</span>;<br>
alpha=mexProximalFlat(U,param);<br>
<br>
<span class="c001">% test l1linf</span><br>
<span class="c002">fprintf</span>(<span class="c003">'\nprox mixed norm l1/linf\n'</span>);<br>
param.regul=<span class="c003">'l1linf'</span>;<br>
alpha=mexProximalFlat(U,param);<br>
<br>
<span class="c001">% test l1l2+l1</span><br>
<span class="c002">fprintf</span>(<span class="c003">'\nprox mixed norm l1/l2 + l1\n'</span>);<br>
param.regul=<span class="c003">'l1l2+l1'</span>;<br>
param.lambda2=0.1;<br>
alpha=mexProximalFlat(U,param);<br>
<br>
<span class="c001">% test l1linf+l1</span><br>
<span class="c002">fprintf</span>(<span class="c003">'\nprox mixed norm l1/linf + l1\n'</span>);<br>
param.regul=<span class="c003">'l1linf+l1'</span>;<br>
param.lambda2=0.1;<br>
alpha=mexProximalFlat(U,param);<br>
<br>
<span class="c001">% test l1linf-row-column</span><br>
<span class="c002">fprintf</span>(<span class="c003">'\nprox mixed norm l1/linf on rows and columns\n'</span>);<br>
param.regul=<span class="c003">'l1linf-row-column'</span>;<br>
param.lambda2=0.1;<br>
alpha=mexProximalFlat(U,param);<br>
<br>
<span class="c001">% test none</span><br>
<span class="c002">fprintf</span>(<span class="c003">'\nprox no regularization\n'</span>);<br>
param.regul=<span class="c003">'none'</span>;<br>
alpha=mexProximalFlat(U,param);</td></tr>
</table>
<h3 class="subsection" id="sec26">5.3  Function mexProximalTree</h3>
<p>
This function computes the proximal operators associated to tree-structured regularization functions, for input signals <span class="c009">U</span>=[<span class="c009">u</span><sup>1</sup>,…,<span class="c009">u</span><sup><span class="c007">n</span></sup>] in ℝ<sup><span class="c007">p</span> × <span class="c007">n</span></sup>, and a tree-structured set of groups [<a href="doc_spams010.html#jenatton3">15</a>], it computes a matrix <span class="c009">V</span>=[<span class="c009">v</span><sup>1</sup>,…,<span class="c009">v</span><sup><span class="c007">n</span></sup>] in ℝ<sup><span class="c007">p</span> × <span class="c007">n</span></sup>. When one uses a regularization function on vectors, it computes a column <span class="c009">v</span> of <span class="c009">V</span> for every column <span class="c009">u</span> of <span class="c009">U</span>:
</p><table class="display dcenter"><tr class="c016"><td class="dcell">
</td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">v</span> ∈ ℝ<sup><span class="c007">p</span></sup></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012">2</td></tr>
</table></td><td class="dcell">||<span class="c009">u</span>−<span class="c009">v</span>||<sub>2</sub><sup>2</sup> + λ </td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">g</span> ∈ <span class="c008">T</span></td></tr>
</table></td><td class="dcell"> η<sup><span class="c007">g</span></sup> ||<span class="c009">v</span><sub><span class="c007">g</span></sub>||<sub>2</sub>,
    (40)</td></tr>
</table><p>
or
</p><table class="display dcenter"><tr class="c016"><td class="dcell">
</td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">v</span> ∈ ℝ<sup><span class="c007">p</span></sup></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012">2</td></tr>
</table></td><td class="dcell">||<span class="c009">u</span>−<span class="c009">v</span>||<sub>2</sub><sup>2</sup> + λ </td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">g</span> ∈ <span class="c008">T</span></td></tr>
</table></td><td class="dcell"> η<sup><span class="c007">g</span></sup> ||<span class="c009">v</span><sub><span class="c007">g</span></sub>||<sub>∞</sub>,
    (41)</td></tr>
</table><p>
or
</p><table class="display dcenter"><tr class="c016"><td class="dcell">
</td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">v</span> ∈ ℝ<sup><span class="c007">p</span></sup></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012">2</td></tr>
</table></td><td class="dcell">||<span class="c009">u</span>−<span class="c009">v</span>||<sub>2</sub><sup>2</sup> + λ </td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">g</span> ∈ <span class="c008">T</span></td></tr>
</table></td><td class="dcell"> δ<sup><span class="c007">g</span></sup>(<span class="c009">v</span>),
    (42)</td></tr>
</table><p>
where δ<sup><span class="c007">g</span></sup>(<span class="c009">v</span>)=0 if <span class="c009">v</span><sub><span class="c007">g</span></sub>=0 and 1 otherwise (see appendix of [<a href="doc_spams010.html#jenatton4">16</a>]).</p><p>When the multi-task tree-structured regularization function is used, it solves
</p><table class="display dcenter"><tr class="c016"><td class="dcell">
</td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">V</span> ∈ ℝ<sup><span class="c007">p</span>× <span class="c007">n</span></sup></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012">2</td></tr>
</table></td><td class="dcell">||<span class="c009">U</span>−<span class="c009">V</span>||<sub><span class="c007">F</span></sub><sup>2</sup> + λ </td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c007">n</span></td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">i</span>=1</td></tr>
</table></td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">g</span> ∈ <span class="c008">T</span></td></tr>
</table></td><td class="dcell"> η<sup><span class="c007">g</span></sup> ||<span class="c009">v</span><sub><span class="c007">g</span></sub><sup><span class="c007">i</span></sup>||<sub>∞</sub>+ λ<sub>2</sub></td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">g</span> ∈ <span class="c008">T</span></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">max</td></tr>
<tr><td class="dcell c012"><span class="c007">j</span> ∈ <span class="c007">g</span></td></tr>
</table></td><td class="dcell"> ||<span class="c009">v</span><sub><span class="c007">g</span></sub><sup><span class="c007">j</span></sup>||<sub>∞</sub>,
    (43)</td></tr>
</table><p> 
which is a formulation presented in [<a href="doc_spams010.html#mairal10">22</a>].</p><p>This function can also be used for computing the proximal operators addressed by mexProximalFlat (it will just not take into account the tree structure). The way the tree is incoded is presented below, (and examples are given in the file test_ProximalTree.m, with more usage details:
</p><table class="lstframe c011"><tr><td class="mouselstlisting"><span class="c001">% <br>
% Usage:  [V [val_regularizer]]=mexProximalTree(U,tree,param);<br>
%<br>
% Name: mexProximalTree<br>
%<br>
% Description: mexProximalTree computes a proximal operator. Depending<br>
%         on the value of param.regul, it computes <br>
%<br>
%         Given an input matrix U=[u^1,\ldots,u^n], and a tree-structured set of groups T,<br>
%         it returns a matrix V=[v^1,\ldots,v^n]:<br>
%         <br>
%         when the regularization function is for vectors,<br>
%         for every column u of U, it compute a column v of V solving<br>
%         if param.regul='tree-l0'<br>
%             argmin 0.5||u-v||_2^2 + lambda \sum_{g \in T} \delta^g(v)<br>
%         if param.regul='tree-l2'<br>
%           for all i, v^i = <br>
%             argmin 0.5||u-v||_2^2 + lambda\sum_{g \in T} \eta_g||v_g||_2<br>
%         if param.regul='tree-linf'<br>
%           for all i, v^i = <br>
%             argmin 0.5||u-v||_2^2 + lambda\sum_{g \in T} \eta_g||v_g||_inf<br>
%<br>
%         when the regularization function is for matrices:<br>
%         if param.regul='multi-task-tree'<br>
%            V=argmin 0.5||U-V||_F^2 + lambda \sum_{i=1}^n\sum_{g \in T} \eta_g||v^i_g||_inf + ...<br>
%                                                lambda_2 \sum_{g \in T} \eta_g max_{j in g}||V_j||_{inf}<br>
%         <br>
%         it can also be used with any non-tree-structured regularization addressed by mexProximalFlat<br>
%<br>
%         for all these regularizations, it is possible to enforce non-negativity constraints<br>
%         with the option param.pos, and to prevent the last row of U to be regularized, with<br>
%         the option param.intercept<br>
%<br>
% Inputs: U:  double m x n matrix   (input signals)<br>
%               m is the signal size<br>
%         tree: struct<br>
%               with four fields, eta_g, groups, own_variables and N_own_variables.<br>
%<br>
%               The tree structure requires a particular organization of groups and variables<br>
%                  * Let us denote by N = |T|, the number of groups.<br>
%                    the groups should be ordered T={g1,g2,\ldots,gN} such that if gi is included<br>
%                    in gj, then j &lt;= i. g1 should be the group at the root of the tree <br>
%                    and contains every variable.<br>
%                  * Every group is a set of  contiguous indices for instance <br>
%                    gi={3,4,5} or gi={4,5,6,7} or gi={4}, but not {3,5};<br>
%                  * We define root(gi) as the indices of the variables that are in gi,<br>
%                    but not in its descendants. For instance for<br>
%                    T={ g1={1,2,3,4},g2={2,3},g3={4} }, then, root(g1)={1}, <br>
%                    root(g2)={2,3}, root(g3)={4},<br>
%                    We assume that for all i, root(gi) is a set of contigous variables<br>
%                  * We assume that the smallest of root(gi) is also the smallest index of gi.<br>
%<br>
%                  For instance, <br>
%                    T={ g1={1,2,3,4},g2={2,3},g3={4} }, is a valid set of groups.<br>
%                    but we can not have<br>
%                    T={ g1={1,2,3,4},g2={1,2},g3={3} }, since root(g1)={4} and 4 is not the<br>
%                    smallest element in g1.<br>
%<br>
%               We do not lose generality with these assumptions since they can be fullfilled for any<br>
%               tree-structured set of groups after a permutation of variables and a correct ordering of the<br>
%               groups.<br>
%               see more examples in test_ProximalTree.m of valid tree-structured sets of groups.<br>
%               <br>
%               The first fields sets the weights for every group<br>
%                  tree.eta_g            double N vector <br>
%  <br>
%               The next field sets inclusion relations between groups <br>
%               (but not between groups and variables):<br>
%                  tree.groups           sparse (double or boolean) N x N matrix  <br>
%                  the (i,j) entry is non-zero if and only if i is different than j and <br>
%                  gi is included in gj.<br>
%                  the first column corresponds to the group at the root of the tree.<br>
%<br>
%               The next field define the smallest index of each group gi, <br>
%               which is also the smallest index of root(gi)<br>
%               tree.own_variables    int32 N vector<br>
%<br>
%               The next field define for each group gi, the size of root(gi)<br>
%               tree.N_own_variables  int32 N vector <br>
%<br>
%               examples are given in test_ProximalTree.m<br>
%<br>
%         param: struct<br>
%               param.lambda  (regularization parameter)<br>
%               param.regul (choice of regularization, see above)<br>
%               param.lambda2  (optional, regularization parameter)<br>
%               param.lambda3  (optional, regularization parameter)<br>
%               param.verbose (optional, verbosity level, false by default)<br>
%               param.intercept (optional, last row of U is not regularized,<br>
%                 false by default)<br>
%               param.pos (optional, adds positivity constraints on the<br>
%                 coefficients, false by default)<br>
%               param.transpose (optional, transpose the matrix in the regularization function)<br>
%               param.size_group (optional, for regularization functions assuming a group<br>
%                 structure). It is a scalar. When param.groups is not specified, it assumes<br>
%                 that the groups are the sets of consecutive elements of size param.size_group<br>
%               param.numThreads (optional, number of threads for exploiting<br>
%                 multi-core / multi-cpus. By default, it takes the value -1,<br>
%                 which automatically selects all the available CPUs/cores).<br>
%<br>
% Output: V: double m x n matrix (output coefficients)<br>
%         val_regularizer: double 1 x n vector (value of the regularization<br>
%         term at the optimum).<br>
%<br>
%<br>
% Author: Julien Mairal, 2010</span></td></tr>
</table><p>
The following piece of code illustrates how to use this function.
</p><table class="lstframe c011"><tr><td class="mouselstlisting">U=<span class="c002">randn</span>(10,1000);<br>
<br>
param.lambda=0.1; <span class="c001">% regularization parameter</span><br>
param.num_threads=-1; <span class="c001">% all cores (-1 by default)</span><br>
param.verbose=true;   <span class="c001">% verbosity, false by default</span><br>
param.pos=false;       <span class="c001">% can be used with all the other regularizations</span><br>
param.intercept=false; <span class="c001">% can be used with all the other regularizations     </span><br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'First tree example\n'</span>);<br>
<span class="c001">% Example 1 of tree structure<br>
% tree structured groups:<br>
% g1= {0 1 2 3 4 5 6 7 8 9}<br>
% g2= {2 3 4}<br>
% g3= {5 6 7 8 9}</span><br>
tree.own_variables=int32([0 2 5]);   <span class="c001">% pointer to the first variable of each group</span><br>
tree.N_own_variables=int32([2 3 5]); <span class="c001">% number of "root" variables in each group</span><br>
                              <span class="c001">% (variables that are in a group, but not in its descendants).</span><br>
                              <span class="c001">% for instance root(g1)={0,1}, root(g2)={2 3 4}, root(g3)={5 6 7 8 9}</span><br>
tree.eta_g=[1 1 1];           <span class="c001">% weights for each group, they should be non-zero to use fenchel duality</span><br>
tree.groups=<span class="c002">sparse</span>([0 0 0; ...<br>
                    1 0 0; ...<br>
                    1 0 0]);    <span class="c001">% first group should always be the root of the tree</span><br>
                                <span class="c001">% non-zero entriees mean inclusion relation ship, here g2 is a children of g1,</span><br>
                                <span class="c001">% g3 is a children of g1</span><br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\ntest prox tree-l0\n'</span>);<br>
param.regul=<span class="c003">'tree-l0'</span>;<br>
alpha=mexProximalTree(U,tree,param);<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\ntest prox tree-l2\n'</span>);<br>
param.regul=<span class="c003">'tree-l2'</span>;<br>
alpha=mexProximalTree(U,tree,param);<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\ntest prox tree-linf\n'</span>);<br>
param.regul=<span class="c003">'tree-linf'</span>;<br>
alpha=mexProximalTree(U,tree,param);<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'Second tree example\n'</span>);<br>
<span class="c001">% Example 2 of tree structure<br>
% tree structured groups:<br>
% g1= {0 1 2 3 4 5 6 7 8 9}    root(g1) = { };<br>
% g2= {0 1 2 3 4 5}            root(g2) = {0 1 2};<br>
% g3= {3 4}                    root(g3) = {3 4};<br>
% g4= {5}                      root(g4) = {5};<br>
% g5= {6 7 8 9}                root(g5) = { };<br>
% g6= {6 7}                    root(g6) = {6 7};<br>
% g7= {8 9}                    root(g7) = {8};<br>
% g8 = {9}                     root(g8) = {9};</span><br>
tree.own_variables=  int32([0 0 3 5 6 6 8 9]);   <span class="c001">% pointer to the first variable of each group</span><br>
tree.N_own_variables=int32([0 3 2 1 0 2 1 1]); <span class="c001">% number of "root" variables in each group</span><br>
tree.eta_g=[1 1 1 2 2 2 2.5 2.5];<br>
tree.groups=<span class="c002">sparse</span>([0 0 0 0 0 0 0 0; ...<br>
                    1 0 0 0 0 0 0 0; ...<br>
                    0 1 0 0 0 0 0 0; ...<br>
                    0 1 0 0 0 0 0 0; ...<br>
                    1 0 0 0 0 0 0 0; ...<br>
                    0 0 0 0 1 0 0 0; ...<br>
                    0 0 0 0 1 0 0 0; ...<br>
                    0 0 0 0 0 0 1 0]);  <span class="c001">% first group should always be the root of the tree</span><br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\ntest prox tree-l0\n'</span>);<br>
param.regul=<span class="c003">'tree-l0'</span>;<br>
alpha=mexProximalTree(U,tree,param);<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\ntest prox tree-l2\n'</span>);<br>
param.regul=<span class="c003">'tree-l2'</span>;<br>
alpha=mexProximalTree(U,tree,param);<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\ntest prox tree-linf\n'</span>);<br>
param.regul=<span class="c003">'tree-linf'</span>;<br>
alpha=mexProximalTree(U,tree,param);<br>
<br>
<span class="c001">% mexProximalTree also works with non-tree-structured regularization functions</span><br>
<span class="c002">fprintf</span>(<span class="c003">'\nprox l1, intercept, positivity constraint\n'</span>);<br>
param.regul=<span class="c003">'l1'</span>;<br>
param.pos=true;       <span class="c001">% can be used with all the other regularizations</span><br>
param.intercept=true; <span class="c001">% can be used with all the other regularizations     </span><br>
alpha=mexProximalTree([U; ones(1,<span class="c002">size</span>(U,2))],tree,param);<br>
<br>
<span class="c001">% Example of multi-task tree</span><br>
<span class="c002">fprintf</span>(<span class="c003">'\nprox multi-task tree\n'</span>);<br>
param.pos=false;<br>
param.intercept=false;<br>
param.lambda2=param.lambda;<br>
param.regul=<span class="c003">'multi-task-tree'</span>;  <span class="c001">% with linf</span><br>
alpha=mexProximalTree(U,tree,param);<br>
<br>
<br>
tree.own_variables=int32([0 1 2 3 4 5 6]);   <span class="c001">% pointer to the first variable of each group</span><br>
tree.N_own_variables=int32([1 1 1 1 1 1]); <span class="c001">% number of "root" variables in each group</span><br>
                              <span class="c001">% (variables that are in a group, but not in its descendants).</span><br>
                              <span class="c001">% for instance root(g1)={0,1}, root(g2)={2 3 4}, root(g3)={5 6 7 8 9}</span><br>
tree.eta_g=[1 1 1 1 1 1];           <span class="c001">% weights for each group, they should be non-zero to use fenchel duality</span><br>
tree.groups=<span class="c002">sparse</span>([0 0 0; ...<br>
                    1 0 0; ...<br>
                    1 0 0]);    <span class="c001">% first group should always be the root of the tree</span><br>
                                <span class="c001">% non-zero entriees mean inclusion relation ship, here g2 is a children of g1,</span><br>
                                <span class="c001">% g3 is a children of g1</span></td></tr>
</table>
<h3 class="subsection" id="sec27">5.4  Function mexProximalGraph</h3>
<p>
This function computes the proximal operators associated to structured sparse regularization, for input signals <span class="c009">U</span>=[<span class="c009">u</span><sup>1</sup>,…,<span class="c009">u</span><sup><span class="c007">n</span></sup>] in ℝ<sup><span class="c007">p</span> × <span class="c007">n</span></sup>, and a set of groups [<a href="doc_spams010.html#mairal10">22</a>], it returns a matrix <span class="c009">V</span>=[<span class="c009">v</span><sup>1</sup>,…,<span class="c009">v</span><sup><span class="c007">n</span></sup>] in ℝ<sup><span class="c007">p</span> × <span class="c007">n</span></sup>.
When one uses a regularization function on vectors, it computes a column <span class="c009">v</span> of <span class="c009">V</span> for every column <span class="c009">u</span> of <span class="c009">U</span>:
</p><table class="display dcenter"><tr class="c016"><td class="dcell">
</td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">v</span> ∈ ℝ<sup><span class="c007">p</span></sup></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012">2</td></tr>
</table></td><td class="dcell">||<span class="c009">u</span>−<span class="c009">v</span>||<sub>2</sub><sup>2</sup> + λ </td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">g</span> ∈ <span class="c008">G</span></td></tr>
</table></td><td class="dcell"> η<sup><span class="c007">g</span></sup> ||<span class="c009">v</span><sub><span class="c007">g</span></sub>||<sub>∞</sub>,
    (44)</td></tr>
</table><p>
or with a regularization function on matrices, it computes <span class="c009">V</span> solving
</p><table class="display dcenter"><tr class="c016"><td class="dcell">
</td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">V</span> ∈ ℝ<sup><span class="c007">p</span>× <span class="c007">n</span></sup></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012">2</td></tr>
</table></td><td class="dcell">||<span class="c009">U</span>−<span class="c009">V</span>||<sub><span class="c007">F</span></sub><sup>2</sup> + λ </td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c007">n</span></td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">i</span>=1</td></tr>
</table></td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">g</span> ∈ <span class="c008">G</span></td></tr>
</table></td><td class="dcell"> η<sup><span class="c007">g</span></sup> ||<span class="c009">v</span><sub><span class="c007">g</span></sub><sup><span class="c007">i</span></sup>||<sub>∞</sub>+ λ<sub>2</sub></td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">g</span> ∈ <span class="c008">G</span></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">max</td></tr>
<tr><td class="dcell c012"><span class="c007">j</span> ∈ <span class="c007">g</span></td></tr>
</table></td><td class="dcell"> ||<span class="c009">v</span><sub><span class="c007">g</span></sub><sup><span class="c007">j</span></sup>||<sub>∞</sub>,
    (45)</td></tr>
</table><p> 
This function can also be used for computing the proximal operators addressed by mexProximalFlat. The way the graph is incoded is presented below (and also in the example file test_ProximalGraph.m, with more usage details:
</p><table class="lstframe c011"><tr><td class="mouselstlisting"><span class="c001">% <br>
% Usage:   [V [val_regularizer]]=mexProximalGraph(U,graph,param);<br>
%<br>
% Name: mexProximalGraph<br>
%<br>
% Description: mexProximalGraph computes a proximal operator. Depending<br>
%         on the value of param.regul, it computes <br>
%<br>
%         Given an input matrix U=[u^1,\ldots,u^n], and a set of groups G,<br>
%         it computes a matrix V=[v^1,\ldots,v^n] such that<br>
%<br>
%         if param.regul='graph'<br>
%         for every column u of U, it computes a column v of V solving<br>
%             argmin 0.5||u-v||_2^2 + lambda\sum_{g \in G} \eta_g||v_g||_inf<br>
%<br>
%         if param.regul='graph+ridge'<br>
%         for every column u of U, it computes a column v of V solving<br>
%             argmin 0.5||u-v||_2^2 + lambda\sum_{g \in G} \eta_g||v_g||_inf + lambda_2||v||_2^2<br>
%<br>
%<br>
%         if param.regul='multi-task-graph'<br>
%            V=argmin 0.5||U-V||_F^2 + lambda \sum_{i=1}^n\sum_{g \in G} \eta_g||v^i_g||_inf + ...<br>
%                                                lambda_2 \sum_{g \in G} \eta_g max_{j in g}||V_j||_{inf}<br>
%         <br>
%         it can also be used with any regularization addressed by mexProximalFlat<br>
%<br>
%         for all these regularizations, it is possible to enforce non-negativity constraints<br>
%         with the option param.pos, and to prevent the last row of U to be regularized, with<br>
%         the option param.intercept<br>
%<br>
% Inputs: U:  double p x n matrix   (input signals)<br>
%               m is the signal size<br>
%         graph: struct<br>
%               with three fields, eta_g, groups, and groups_var<br>
%<br>
%               The first fields sets the weights for every group<br>
%                  graph.eta_g            double N vector <br>
%  <br>
%               The next field sets inclusion relations between groups <br>
%               (but not between groups and variables):<br>
%                  graph.groups           sparse (double or boolean) N x N matrix  <br>
%                  the (i,j) entry is non-zero if and only if i is different than j and <br>
%                  gi is included in gj.<br>
%               <br>
%               The next field sets inclusion relations between groups and variables<br>
%                  graph.groups_var       sparse (double or boolean) p x N matrix<br>
%                  the (i,j) entry is non-zero if and only if the variable i is included <br>
%                  in gj, but not in any children of gj.<br>
%<br>
%               examples are given in test_ProximalGraph.m<br>
%<br>
%         param: struct<br>
%               param.lambda  (regularization parameter)<br>
%               param.regul (choice of regularization, see above)<br>
%               param.lambda2  (optional, regularization parameter)<br>
%               param.lambda3  (optional, regularization parameter)<br>
%               param.verbose (optional, verbosity level, false by default)<br>
%               param.intercept (optional, last row of U is not regularized,<br>
%                 false by default)<br>
%               param.pos (optional, adds positivity constraints on the<br>
%                 coefficients, false by default)<br>
%               param.numThreads (optional, number of threads for exploiting<br>
%                 multi-core / multi-cpus. By default, it takes the value -1,<br>
%                 which automatically selects all the available CPUs/cores).<br>
%<br>
% Output: V: double p x n matrix (output coefficients)<br>
%         val_regularizer: double 1 x n vector (value of the regularization<br>
%         term at the optimum).<br>
%<br>
% Author: Julien Mairal, 2010</span></td></tr>
</table><p>
The following piece of code illustrates how to use this function.
</p><table class="lstframe c011"><tr><td class="mouselstlisting">U=<span class="c002">randn</span>(10,1000);<br>
<br>
param.lambda=0.1; <span class="c001">% regularization parameter</span><br>
param.num_threads=-1; <span class="c001">% all cores (-1 by default)</span><br>
param.verbose=true;   <span class="c001">% verbosity, false by default</span><br>
param.pos=false;       <span class="c001">% can be used with all the other regularizations</span><br>
param.intercept=false; <span class="c001">% can be used with all the other regularizations     </span><br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'First graph example\n'</span>);<br>
<span class="c001">% Example 1 of graph structure<br>
% groups:<br>
% g1= {0 1 2 3}<br>
% g2= {3 4 5 6}<br>
% g3= {6 7 8 9}</span><br>
graph.eta_g=[1 1 1];<br>
graph.groups=<span class="c002">sparse</span>(<span class="c002">zeros</span>(3));<br>
graph.groups_var=<span class="c002">sparse</span>([1 0 0;<br>
                         1 0 0;<br>
                         1 0 0;<br>
                         1 1 0;<br>
                         0 1 0;<br>
                         0 1 0;<br>
                         0 1 1;<br>
                         0 0 1;<br>
                         0 0 1;<br>
                         0 0 1]);<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\ntest prox graph\n'</span>);<br>
param.regul=<span class="c003">'graph'</span>;<br>
alpha=mexProximalGraph(U,graph,param);<br>
<br>
<span class="c001">% Example 2 of graph structure<br>
% groups:<br>
% g1= {0 1 2 3}<br>
% g2= {3 4 5 6}<br>
% g3= {6 7 8 9}<br>
% g4= {0 1 2 3 4 5}<br>
% g5= {6 7 8}</span><br>
graph.eta_g=[1 1 1 1 1];<br>
graph.groups=<span class="c002">sparse</span>([0 0 0 1 0;<br>
                     0 0 0 0 0;<br>
                     0 0 0 0 0;<br>
                     0 0 0 0 0;<br>
                     0 0 1 0 0]);   <span class="c001">% g5 is included in g3, and g2 is included in g4</span><br>
graph.groups_var=<span class="c002">sparse</span>([1 0 0 0 0;<br>
                         1 0 0 0 0;<br>
                         1 0 0 0 0 ;<br>
                         1 1 0 0 0;<br>
                         0 1 0 1 0;<br>
                         0 1 0 1 0;<br>
                         0 1 0 0 1;<br>
                         0 0 0 0 1;<br>
                         0 0 0 0 1;<br>
                         0 0 1 0 0]); <span class="c001">% represents direct inclusion relations </span><br>
                                      <span class="c001">% between groups (columns) and variables (rows)</span><br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\ntest prox graph\n'</span>);<br>
param.regul=<span class="c003">'graph'</span>;<br>
alpha=mexProximalGraph(U,graph,param);<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\ntest prox multi-task-graph\n'</span>);<br>
param.regul=<span class="c003">'multi-task-graph'</span>;<br>
param.lambda2=0.1;<br>
alpha=mexProximalGraph(U,graph,param);<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\ntest no regularization\n'</span>);<br>
param.regul=<span class="c003">'none'</span>;<br>
alpha=mexProximalGraph(U,graph,param);</td></tr>
</table>
<h3 class="subsection" id="sec28">5.5  Function mexProximalPathCoding</h3>
<p>
This function computes the proximal operators associated to the path coding penalties of [<a href="doc_spams010.html#mairal14">24</a>]. 
</p><table class="lstframe c011"><tr><td class="mouselstlisting"><span class="c001">% <br>
% Usage:   [V [val_regularizer]]=mexProximalPathCoding(U,DAG,param);<br>
%<br>
% Name: mexProximalPathCoding<br>
%<br>
% Description: mexProximalPathCoding computes a proximal operator for<br>
%         the path coding penalties of http://arxiv.org/abs/1204.4539<br>
%<br>
%         Given an input matrix U=[u^1,\ldots,u^n], <br>
%<br>
%<br>
% Inputs: U:  double p x n matrix   (input signals)<br>
%               m is the signal size<br>
%         DAG:  struct<br>
%               with three fields, weights, start_weights, stop_weights<br>
%         for a graph with |V| nodes and |E| arcs,<br>
%         DAG.weights: sparse double |V| x |V| matrix. Adjacency<br>
%               matrix. The non-zero entries represent costs on arcs<br>
%               linking two nodes.<br>
%         DAG.start_weights: dense double |V| vector. Represent the costs<br>
%               of starting a path from a specific node.<br>
%         DAG.stop_weights: dense double |V| vector. Represent the costs<br>
%               of ending a path at a specific node.<br>
%<br>
%         if param.regul='graph-path-l0', non-convex penalty<br>
%         if param.regul='graph-path-conv', convex penalty<br>
%<br>
%         param: struct<br>
%               param.lambda  (regularization parameter)<br>
%               param.regul (choice of regularization, see above)<br>
%               param.verbose (optional, verbosity level, false by default)<br>
%               param.intercept (optional, last row of U is not regularized,<br>
%                 false by default)<br>
%               param.pos (optional, adds positivity constraints on the<br>
%                 coefficients, false by default)<br>
%               param.precision (optional, by default a very large integer.<br>
%                 It returns approximate proximal operator by choosing a small integer,<br>
%                 for example, 100 or 1000.<br>
%               param.numThreads (optional, number of threads for exploiting<br>
%                 multi-core / multi-cpus. By default, it takes the value -1,<br>
%                 which automatically selects all the available CPUs/cores).<br>
%<br>
% Output: V: double p x n matrix (output coefficients)<br>
%         val_regularizer: double 1 x n vector (value of the regularization<br>
%         term at the optimum).<br>
%<br>
% Author: Julien Mairal, 2012</span></td></tr>
</table><p>
The following piece of code illustrates how to use this function.
</p><table class="lstframe c011"><tr><td class="mouselstlisting"><span class="c002">clear</span> <span class="c002">all</span>;<br>
<span class="c002">rand</span>(<span class="c003">'seed'</span>,0);<br>
<span class="c002">randn</span>(<span class="c003">'seed'</span>,0);<br>
<span class="c002">fprintf</span>(<span class="c003">'test mexProximalPathCoding\n'</span>);<br>
p=100;<br>
<span class="c001">% generate a DAG</span><br>
G=sprand(p,p,0.02);<br>
G=mexRemoveCyclesGraph(G);<br>
<span class="c002">fprintf</span>(<span class="c003">'\n'</span>);<br>
<br>
<span class="c001">% generate a data matrix</span><br>
U=<span class="c002">randn</span>(p,10);<br>
U=U-<span class="c002">mean</span>(U(:));<br>
U=mexNormalize(U);<br>
<br>
<span class="c001">% input graph</span><br>
graph.weights=G;<br>
graph.stop_weights=<span class="c002">zeros</span>(1,p);<br>
graph.start_weights=10*ones(1,p);<br>
<br>
<span class="c001">% FISTA parameters</span><br>
param.num_threads=-1; <span class="c001">% all cores (-1 by default)</span><br>
param.verbose=true;   <span class="c001">% verbosity, false by default</span><br>
param.lambda=0.05; <span class="c001">% regularization parameter</span><br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'Proximal convex path penalty\n'</span>);<br>
param.regul=<span class="c003">'graph-path-conv'</span>;<br>
<span class="c002">tic</span><br>
[V1 optim]=mexProximalPathCoding(U,graph,param);<br>
t=<span class="c002">toc</span>;<br>
num=mexCountConnexComponents(graph.weights,V1(:,1));<br>
<span class="c002">fprintf</span>(<span class="c003">'Num of connected components: %d\n'</span>,num);<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'Proximal non-convex path penalty\n'</span>);<br>
param.regul=<span class="c003">'graph-path-l0'</span>;<br>
param.lambda=0.005;<br>
<span class="c002">tic</span><br>
[V2 optim]=mexProximalPathCoding(U,graph,param);<br>
t=<span class="c002">toc</span>;<br>
num=mexCountConnexComponents(graph.weights,V2(:,1));<br>
<span class="c002">fprintf</span>(<span class="c003">'Num of connected components: %d\n'</span>,num);<br>
<br>
graph.start_weights=1*ones(1,p);<br>
param.lambda=0.05;<br>
<span class="c002">fprintf</span>(<span class="c003">'Proximal convex path penalty\n'</span>);<br>
param.regul=<span class="c003">'graph-path-conv'</span>;<br>
<span class="c002">tic</span><br>
[V1 optim]=mexProximalPathCoding(U,graph,param);<br>
t=<span class="c002">toc</span>;<br>
num=mexCountConnexComponents(graph.weights,V1(:,1));<br>
<span class="c002">fprintf</span>(<span class="c003">'Num of connected components: %d\n'</span>,num);<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'Proximal non-convex path penalty\n'</span>);<br>
param.regul=<span class="c003">'graph-path-l0'</span>;<br>
param.lambda=0.005;<br>
<span class="c002">tic</span><br>
[V2 optim]=mexProximalPathCoding(U,graph,param);<br>
t=<span class="c002">toc</span>;<br>
num=mexCountConnexComponents(graph.weights,V2(:,1));<br>
<span class="c002">fprintf</span>(<span class="c003">'Num of connected components: %d\n'</span>,num);</td></tr>
</table><p>
This function is associated to a function to evaluate the penalties:
</p>
<h3 class="subsection" id="sec29">5.6  Function mexEvalPathCoding</h3>
<table class="lstframe c011"><tr><td class="mouselstlisting"><span class="c001">% <br>
% Usage:   [val [paths]]=mexEvalPathCoding(U,DAG,param);<br>
%<br>
% Name: mexEvalPathCoding<br>
%<br>
% Description: mexEvalPathCoding evaluate the path coding penalies <br>
%         of http://arxiv.org/abs/1204.4539 and provides a path <br>
%         decomposition of a vector W.<br>
%<br>
%         Given an input matrix U=[u^1,\ldots,u^n], <br>
%<br>
%<br>
% Inputs: U:  double p x n matrix   (input signals)<br>
%               m is the signal size<br>
%         DAG:  struct<br>
%               with three fields, weights, start_weights, stop_weights<br>
%         for a graph with |V| nodes and |E| arcs,<br>
%         DAG.weights: sparse double |V| x |V| matrix. Adjacency<br>
%               matrix. The non-zero entries represent costs on arcs<br>
%               linking two nodes.<br>
%         DAG.start_weights: dense double |V| vector. Represent the costs<br>
%               of starting a path from a specific node.<br>
%         DAG.stop_weights: dense double |V| vector. Represent the costs<br>
%               of ending a path at a specific node.<br>
%<br>
%         if param.regul='graph-path-l0', non-convex penalty<br>
%         if param.regul='graph-path-conv', convex penalty<br>
%<br>
%         param: struct<br>
%               param.regul (choice of regularization, see above)<br>
%               param.verbose (optional, verbosity level, false by default)<br>
%               param.precision (optional, by default a very large integer.<br>
%                 It returns approximate proximal operator by choosing a small integer,<br>
%                 for example, 100 or 1000.<br>
%               param.numThreads (optional, number of threads for exploiting<br>
%                 multi-core / multi-cpus. By default, it takes the value -1,<br>
%                 which automatically selects all the available CPUs/cores).<br>
%<br>
% Output: V: double 1 x n vector (values of the objective function)<br>
%         paths: optional, double sparse p x k matrix. selected paths for the <br>
%                first column of U<br>
%<br>
% Author: Julien Mairal, 2012</span></td></tr>
</table><p>
The following piece of code illustrates how to use this function.
</p><table class="lstframe c011"><tr><td class="mouselstlisting"><span class="c002">clear</span> <span class="c002">all</span>;<br>
<span class="c002">rand</span>(<span class="c003">'seed'</span>,0);<br>
<span class="c002">randn</span>(<span class="c003">'seed'</span>,0);<br>
<br>
p=100;<br>
G=sprand(p,p,0.05);<br>
G=mexRemoveCyclesGraph(G);<br>
<span class="c002">fprintf</span>(<span class="c003">'\n'</span>);<br>
<br>
<span class="c001">% input graph</span><br>
graph.weights=G;<br>
graph.stop_weights=<span class="c002">zeros</span>(1,p);<br>
graph.start_weights=10*ones(1,p);<br>
<br>
param.regul=<span class="c003">'graph-path-l0'</span>;<br>
U=<span class="c002">randn</span>(p,10);<br>
U=U-<span class="c002">mean</span>(U(:));<br>
U=mexNormalize(U);<br>
param.lambda=0.005;<br>
[V2 optim]=mexProximalPathCoding(U,graph,param);<br>
[vals paths]=mexEvalPathCoding(U,graph,param);</td></tr>
</table><p>After having presented the regularization terms which our software can handle,
we present the various formulations that we address
</p>
<h3 class="subsection" id="sec30">5.7  Problems Addressed</h3>
<p>
We present here regression or classification formulations and their multi-task variants.
</p>
<h4 class="subsubsection" id="sec31">5.7.1  Regression Problems with the Square Loss</h4>
<p> Given a training set {<span class="c009">x</span><sup><span class="c007">i</span></sup>,<span class="c007">y</span><sub><span class="c007">i</span></sub>}<sub><span class="c007">i</span>=1</sub><sup><span class="c007">n</span></sup>, with <span class="c009">x</span><sup><span class="c007">i</span></sup> ∈ ℝ<sup><span class="c007">p</span></sup> and <span class="c007">y</span><sub><span class="c007">i</span></sub> ∈ ℝ for all <span class="c007">i</span> in [ 1;<span class="c007">n</span> ], we address
</p><table class="display dcenter"><tr class="c016"><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">w</span> ∈ ℝ<sup><span class="c007">p</span></sup>, <span class="c007">b</span> ∈ ℝ</td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c007">n</span></td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">i</span>=1</td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012">2</td></tr>
</table></td><td class="dcell">(<span class="c007">y</span><sub><span class="c007">i</span></sub>−<span class="c009">w</span><sup>⊤</sup><span class="c009">x</span><sup><span class="c007">i</span></sup>−<span class="c007">b</span>)<sup>2</sup> + λψ(<span class="c009">w</span>),
</td></tr>
</table><p>
where <span class="c007">b</span> is an optional variable acting as an “intercept”, which is not regularized, and ψ
can be any of the regularization functions presented above. 
Let us consider the vector <span class="c009">y</span> in ℝ<sup><span class="c007">n</span></sup> that carries the entries <span class="c007">y</span><sub><span class="c007">i</span></sub>. 
The problem without the intercept takes the following form, which we have already
encountered in the previous toolbox, but with different notations:
</p><table class="display dcenter"><tr class="c016"><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">w</span> ∈ ℝ<sup><span class="c007">p</span></sup></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012">2</td></tr>
</table></td><td class="dcell">||<span class="c009">y</span>−<span class="c009">Xw</span>||<sub>2</sub><sup>2</sup> + λψ(<span class="c009">w</span>),
</td></tr>
</table><p>
where the <span class="c009">X</span>=[<span class="c009">x</span><sup><span class="c007">i</span></sup>,…,<span class="c009">x</span><sup><span class="c007">n</span></sup>]<sup><span class="c007">T</span></sup> (the <span class="c009">x</span><sup><span class="c007">i</span></sup>’s are here the rows of <span class="c009">X</span>).
</p>
<h4 class="subsubsection" id="sec32">5.7.2  Classification Problems with the Logistic Loss</h4>
<p>
The next formulation that our software can solve is the regularized logistic regression formulation.
We are again given a training set {<span class="c009">x</span><sup><span class="c007">i</span></sup>,<span class="c007">y</span><sub><span class="c007">i</span></sub>}<sub><span class="c007">i</span>=1</sub><sup><span class="c007">n</span></sup>, with <span class="c009">x</span><sup><span class="c007">i</span></sup> ∈
ℝ<sup><span class="c007">p</span></sup>, but the variables <span class="c007">y</span><sub><span class="c007">i</span></sub> are now in {−1,+1} for all <span class="c007">i</span> in
[ 1;<span class="c007">n</span> ]. The optimization problem we address is
</p><table class="display dcenter"><tr class="c016"><td class="dcell">   </td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">w</span> ∈ ℝ<sup><span class="c007">p</span></sup>, <span class="c007">b</span> ∈ ℝ</td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012"><span class="c007">n</span></td></tr>
</table></td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c007">n</span></td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">i</span>=1</td></tr>
</table></td><td class="dcell"> log(1+<span class="c007">e</span><sup>−<span class="c007">y</span><sub><span class="c007">i</span></sub>(<span class="c009">w</span><sup>⊤</sup><span class="c009">x</span><sup><span class="c007">i</span></sup>+<span class="c007">b</span>)</sup> + λψ(<span class="c009">w</span>),
</td></tr>
</table><p>
with again ψ taken to be one of the regularization function presented above, and <span class="c007">b</span> is an optional intercept.
</p>
<h4 class="subsubsection" id="sec33">5.7.3  Multi-class Classification Problems with the Softmax Loss</h4>
<p>
We have also implemented a multi-class logistic classifier (or softmax).
For a classification problem with <span class="c007">r</span> classes, we are given a training set {<span class="c009">x</span><sup><span class="c007">i</span></sup>,<span class="c007">y</span><sub><span class="c007">i</span></sub>}<sub><span class="c007">i</span>=1</sub><sup><span class="c007">n</span></sup>, where the variables <span class="c009">x</span><sup><span class="c007">i</span></sup> are still vectors in ℝ<sup><span class="c007">p</span></sup>, but the <span class="c007">y</span><sub><span class="c007">i</span></sub>’s have integer values in {1,2,…,<span class="c007">r</span>}. The formulation we address is the following multi-class learning problem
</p><table class="display dcenter"><tr class="c016"><td class="dcell">
</td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">W</span> ∈ ℝ<sup><span class="c007">p</span> × <span class="c007">r</span></sup>, <span class="c009">b</span> ∈ ℝ<sup><span class="c007">r</span></sup></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012"><span class="c007">n</span></td></tr>
</table></td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c007">n</span></td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">i</span>=1</td></tr>
</table></td><td class="dcell"> log</td><td class="dcell">⎛<br>
⎜<br>
⎝</td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c007">r</span></td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">j</span>=1</td></tr>
</table></td><td class="dcell"> <span class="c007">e</span><sup> (<span class="c009">w</span><sup><span class="c007">j</span></sup>−<span class="c009">w</span><sup><span class="c009">y</span><sub><span class="c007">i</span></sub></sup>)<sup>⊤</sup><span class="c009">x</span><sup><span class="c007">i</span></sup> + <span class="c009">b</span><sub><span class="c007">j</span></sub>−<span class="c009">b</span><sub><span class="c009">y</span><sub><span class="c007">i</span></sub></sub></sup></td><td class="dcell">⎞<br>
⎟<br>
⎠</td><td class="dcell">+ λ</td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c007">r</span></td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">j</span>=1</td></tr>
</table></td><td class="dcell">ψ(<span class="c009">w</span><sup><span class="c007">j</span></sup>),<a id="software:eq:class"></a>
    (46)</td></tr>
</table><p>
where <span class="c009">W</span> = [<span class="c009">w</span><sup>1</sup>,…,<span class="c009">w</span><sup><span class="c007">r</span></sup>] and the optional vector <span class="c009">b</span> in ℝ<sup><span class="c007">r</span></sup> carries intercepts for each class.
</p>
<h4 class="subsubsection" id="sec34">5.7.4  Multi-task Regression Problems with the Square Loss</h4>
<p>
We are now considering a problem with <span class="c007">r</span> tasks, and a training set
{<span class="c009">x</span><sup><span class="c007">i</span></sup>,<span class="c009">y</span><sup><span class="c007">i</span></sup>}<sub><span class="c007">i</span>=1</sub><sup><span class="c007">n</span></sup>, where the variables <span class="c009">x</span><sup><span class="c007">i</span></sup> are still vectors in ℝ<sup><span class="c007">p</span></sup>, and <span class="c009">y</span><sup><span class="c007">i</span></sup>
is a vector in ℝ<sup><span class="c007">r</span></sup>. We are looking for <span class="c007">r</span> regression vectors <span class="c009">w</span><sup><span class="c007">j</span></sup>, for <span class="c007">j</span>∈ [ 1;<span class="c007">r</span> ], or equivalently for a matrix <span class="c009">W</span>=[<span class="c009">w</span><sup>1</sup>,…,<span class="c009">w</span><sup><span class="c007">r</span></sup>] in ℝ<sup><span class="c007">p</span> × <span class="c007">r</span></sup>. The formulation we address is the following
multi-task regression problem
</p><table class="display dcenter"><tr class="c016"><td class="dcell">         </td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">W</span> ∈ ℝ<sup><span class="c007">p</span> × <span class="c007">r</span></sup>, <span class="c009">b</span> ∈ ℝ<sup><span class="c007">r</span></sup></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c007">r</span></td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">j</span>=1</td></tr>
</table></td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c007">n</span></td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">i</span>=1</td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012">2</td></tr>
</table></td><td class="dcell">(<span class="c009">y</span><sub><span class="c007">j</span></sub><sup><span class="c007">i</span></sup>−<span class="c009">w</span><sup>⊤</sup><span class="c009">x</span><sup><span class="c007">i</span></sup>−<span class="c009">b</span><sub><span class="c007">j</span></sub>)<sup>2</sup> + λψ(<span class="c009">W</span>),
</td></tr>
</table><p>
where ψ is any of the regularization function on matrices we have presented in the previous section.
Note that by introducing the appropriate variables <span class="c009">Y</span>, the problem without intercept could be equivalently rewritten
</p><table class="display dcenter"><tr class="c016"><td class="dcell">   </td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">W</span> ∈ ℝ<sup><span class="c007">p</span> × <span class="c007">r</span></sup></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012">2</td></tr>
</table></td><td class="dcell"> ||<span class="c009">Y</span>−<span class="c009">XW</span>||<sub>F</sub><sup>2</sup> + λψ(<span class="c009">W</span>).
</td></tr>
</table>
<h4 class="subsubsection" id="sec35">5.7.5  Multi-task Classification Problems with the Logistic Loss</h4>
<p>
The multi-task version of the logistic regression follows the same principle.
We consider <span class="c007">r</span> tasks, and a training set
{<span class="c009">x</span><sup><span class="c007">i</span></sup>,<span class="c009">y</span><sup><span class="c007">i</span></sup>}<sub><span class="c007">i</span>=1</sub><sup><span class="c007">n</span></sup>, with the <span class="c009">x</span><sup><span class="c007">i</span></sup>’s in ℝ<sup><span class="c007">p</span></sup>, and the <span class="c009">y</span><sup><span class="c007">i</span></sup>’s
are vectors in {−1,+1}<sup><span class="c007">r</span></sup>. We look for a matrix <span class="c009">W</span>=[<span class="c009">w</span><sup>1</sup>,…,<span class="c009">w</span><sup><span class="c007">r</span></sup>] in ℝ<sup><span class="c007">p</span> × <span class="c007">r</span></sup>. The formulation is the following
multi-task regression problem
</p><table class="display dcenter"><tr class="c016"><td class="dcell">   </td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">W</span> ∈ ℝ<sup><span class="c007">p</span> × <span class="c007">r</span></sup>, <span class="c009">b</span> ∈ ℝ<sup><span class="c007">r</span></sup></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c007">r</span></td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">j</span>=1</td></tr>
</table></td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012"><span class="c007">n</span></td></tr>
</table></td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c007">n</span></td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">i</span>=1</td></tr>
</table></td><td class="dcell"> log</td><td class="dcell">⎛<br>
⎜<br>
⎝</td><td class="dcell">1+<span class="c007">e</span><sup>−<span class="c009">y</span><sub><span class="c007">j</span></sub><sup><span class="c007">i</span></sup>(<span class="c009">w</span><sup>⊤</sup><span class="c009">x</span><sup><span class="c007">i</span></sup>+<span class="c009">b</span><sub><span class="c007">j</span></sub>)</sup></td><td class="dcell">⎞<br>
⎟<br>
⎠</td><td class="dcell">+ λψ(<span class="c009">W</span>).
</td></tr>
</table>
<h4 class="subsubsection" id="sec36">5.7.6  Multi-task and Multi-class Classification Problems with the Softmax Loss</h4>
<p>
The multi-task/multi-class version directly follows from the formulation of Eq. (<a href="#software%3Aeq%3Aclass">46</a>), but associates with each class a task, and as a consequence, regularizes the matrix <span class="c009">W</span> in a particular way:
</p><table class="display dcenter"><tr class="c016"><td class="dcell">   </td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">W</span> ∈ ℝ<sup><span class="c007">p</span> × <span class="c007">r</span></sup>, <span class="c009">b</span> ∈ ℝ<sup><span class="c007">r</span></sup></td></tr>
</table></td><td class="dcell"> </td><td class="dcell"><table class="display"><tr><td class="dcell c012">1</td></tr>
<tr><td class="hbar"></td></tr>
<tr><td class="dcell c012"><span class="c007">n</span></td></tr>
</table></td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c007">n</span></td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">i</span>=1</td></tr>
</table></td><td class="dcell"> log</td><td class="dcell">⎛<br>
⎜<br>
⎝</td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c007">r</span></td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">j</span>=1</td></tr>
</table></td><td class="dcell"> <span class="c007">e</span><sup> (<span class="c009">w</span><sup><span class="c007">j</span></sup>−<span class="c009">w</span><sup><span class="c009">y</span><sub><span class="c007">i</span></sub></sup>)<sup>⊤</sup><span class="c009">x</span><sup><span class="c007">i</span></sup> + <span class="c009">b</span><sub><span class="c007">j</span></sub>−<span class="c009">b</span><sub><span class="c009">y</span><sub><span class="c007">i</span></sub></sub></sup></td><td class="dcell">⎞<br>
⎟<br>
⎠</td><td class="dcell">+ λψ(<span class="c009">W</span>).
</td></tr>
</table><p>
How duality gaps are computed for any of these formulations is presented in Appendix <a href="doc_spams009.html#appendix">A</a>.
We now present the main functions for solving these problems</p>
<h3 class="subsection" id="sec37">5.8  Function mexFistaFlat</h3>
<p>
Given a matrix <span class="c009">X</span>=[<span class="c009">x</span><sup>1</sup>,…,<span class="c009">x</span><sup><span class="c007">p</span></sup>]<sup><span class="c007">T</span></sup> in ℝ<sup><span class="c007">m</span> × <span class="c007">p</span></sup>, and a matrix <span class="c009">Y</span>=[<span class="c009">y</span><sup>1</sup>,…,<span class="c009">y</span><sup><span class="c007">n</span></sup>], it solves the optimization problems presented in the previous section, with the same regularization functions as mexProximalFlat.
see usage details below:
</p><table class="lstframe c011"><tr><td class="mouselstlisting"><span class="c001">% <br>
% Usage: [W [optim]]=mexFistaFlat(Y,X,W0,param);<br>
%<br>
% Name: mexFistaFlat<br>
%<br>
% Description: mexFistaFlat solves sparse regularized problems.<br>
%         X is a design matrix of size m x p<br>
%         X=[x^1,...,x^n]', where the x_i's are the rows of X<br>
%         Y=[y^1,...,y^n] is a matrix of size m x n<br>
%         It implements the algorithms FISTA, ISTA and subgradient descent.<br>
%         <br>
%           - if param.loss='square' and param.regul is a regularization function for vectors,<br>
%             the entries of Y are real-valued,  W = [w^1,...,w^n] is a matrix of size p x n<br>
%             For all column y of Y, it computes a column w of W such that<br>
%               w = argmin 0.5||y- X w||_2^2 + lambda psi(w)<br>
%<br>
%           - if param.loss='square' and param.regul is a regularization function for matrices<br>
%             the entries of Y are real-valued,  W is a matrix of size p x n. <br>
%             It computes the matrix W such that<br>
%               W = argmin 0.5||Y- X W||_F^2 + lambda psi(W)<br>
%            <br>
%           - param.loss='square-missing' : same as param.loss='square', but handles missing data<br>
%             represented by NaN (not a number) in the matrix Y<br>
%<br>
%           - if param.loss='logistic' and param.regul is a regularization function for vectors,<br>
%             the entries of Y are either -1 or +1, W = [w^1,...,w^n] is a matrix of size p x n<br>
%             For all column y of Y, it computes a column w of W such that<br>
%               w = argmin (1/m)sum_{j=1}^m log(1+e^(-y_j x^j' w)) + lambda psi(w),<br>
%             where x^j is the j-th row of X.<br>
%<br>
%           - if param.loss='logistic' and param.regul is a regularization function for matrices<br>
%             the entries of Y are either -1 or +1, W is a matrix of size p x n<br>
%               W = argmin sum_{i=1}^n(1/m)sum_{j=1}^m log(1+e^(-y^i_j x^j' w^i)) + lambda psi(W)<br>
%<br>
%           - if param.loss='multi-logistic' and param.regul is a regularization function for vectors,<br>
%             the entries of Y are in {0,1,...,N} where N is the total number of classes<br>
%             W = [W^1,...,W^n] is a matrix of size p x Nn, each submatrix W^i is of size p x N<br>
%             for all submatrix WW of W, and column y of Y, it computes<br>
%               WW = argmin (1/m)sum_{j=1}^m log(sum_{j=1}^r e^(x^j'(ww^j-ww^{y_j}))) + lambda sum_{j=1}^N psi(ww^j),<br>
%             where ww^j is the j-th column of WW.<br>
%<br>
%           - if param.loss='multi-logistic' and param.regul is a regularization function for matrices,<br>
%             the entries of Y are in {0,1,...,N} where N is the total number of classes<br>
%             W is a matrix of size p x N, it computes<br>
%               W = argmin (1/m)sum_{j=1}^m log(sum_{j=1}^r e^(x^j'(w^j-w^{y_j}))) + lambda psi(W)<br>
%             where ww^j is the j-th column of WW.<br>
%<br>
%           - param.loss='cur' : useful to perform sparse CUR matrix decompositions, <br>
%               W = argmin 0.5||Y-X*W*X||_F^2 + lambda psi(W)<br>
%<br>
%<br>
%         The function psi are those used by mexProximalFlat (see documentation)<br>
%<br>
%         This function can also handle intercepts (last row of W is not regularized),<br>
%         and/or non-negativity constraints on W, and sparse matrices for X<br>
%<br>
% Inputs: Y:  double dense m x n matrix<br>
%         X:  double dense or sparse m x p matrix   <br>
%         W0:  double dense p x n matrix or p x Nn matrix (for multi-logistic loss)<br>
%              initial guess<br>
%         param: struct<br>
%            param.loss (choice of loss, see above)<br>
%            param.regul (choice of regularization, see function mexProximalFlat)<br>
%            param.lambda (regularization parameter)<br>
%            param.lambda2 (optional, regularization parameter, 0 by default)<br>
%            param.lambda3 (optional, regularization parameter, 0 by default)<br>
%            param.verbose (optional, verbosity level, false by default)<br>
%            param.pos (optional, adds positivity constraints on the<br>
%                coefficients, false by default)<br>
%            param.transpose (optional, transpose the matrix in the regularization function)<br>
%            param.size_group (optional, for regularization functions assuming a group<br>
%                 structure)<br>
%            param.groups (int32, optional, for regularization functions assuming a group<br>
%                 structure, see mexProximalFlat)<br>
%            param.numThreads (optional, number of threads for exploiting<br>
%                multi-core / multi-cpus. By default, it takes the value -1,<br>
%                which automatically selects all the available CPUs/cores).<br>
%            param.max_it (optional, maximum number of iterations, 100 by default)<br>
%            param.it0 (optional, frequency for computing duality gap, every 10 iterations by default)<br>
%            param.tol (optional, tolerance for stopping criteration, which is a relative duality gap<br>
%                if it is available, or a relative change of parameters).<br>
%            param.gamma (optional, multiplier for increasing the parameter L in fista, 1.5 by default)<br>
%            param.L0 (optional, initial parameter L in fista, 0.1 by default, should be small enough)<br>
%            param.fixed_step (deactive the line search for L in fista and use param.L0 instead)<br>
%            param.linesearch_mode (line-search scheme when ista=true:<br>
%                      0: default, monotonic backtracking scheme<br>
%                      1: monotonic backtracking scheme, with restart at each iteration<br>
%                      2: Barzilai-Borwein step sizes (similar to SparSA by Wright et al.)<br>
%                      3: non-monotonic backtracking<br>
%            param.compute_gram (optional, pre-compute X^TX, false by default).<br>
%            param.intercept (optional, do not regularize last row of W, false by default).<br>
%            param.ista (optional, use ista instead of fista, false by default).<br>
%            param.subgrad (optional, if not param.ista, use subradient descent instead of fista, false by default).<br>
%            param.a, param.b (optional, if param.subgrad, the gradient step is a/(t+b)<br>
%            also similar options as mexProximalFlat<br>
%<br>
%            the function also implements the ADMM algorithm via an option param.admm=true. It is not documented<br>
%            and you need to look at the source code to use it.<br>
%<br>
% Output:  W:  double dense p x n matrix or p x Nn matrix (for multi-logistic loss)<br>
%          optim: optional, double dense 4 x n matrix.<br>
%              first row: values of the objective functions.<br>
%              third row: values of the relative duality gap (if available)<br>
%              fourth row: number of iterations<br>
%<br>
% Author: Julien Mairal, 2010</span></td></tr>
</table><p>
The following piece of code illustrates how to use this function.
</p><table class="lstframe c011"><tr><td class="mouselstlisting"><span class="c002">format</span> compact;<br>
<span class="c002">randn</span>(<span class="c003">'seed'</span>,0);<br>
param.numThreads=-1; <span class="c001">% all cores (-1 by default)</span><br>
param.verbose=true;   <span class="c001">% verbosity, false by default</span><br>
param.lambda=0.05; <span class="c001">% regularization parameter</span><br>
param.it0=10;      <span class="c001">% frequency for duality gap computations</span><br>
param.max_it=200; <span class="c001">% maximum number of iterations</span><br>
param.L0=0.1;<br>
param.tol=1e-3;<br>
param.intercept=false;<br>
param.pos=false;<br>
param.ista=false;<br>
<br>
X=<span class="c002">randn</span>(100,200);<br>
X=X-repmat(<span class="c002">mean</span>(X),[<span class="c002">size</span>(X,1) 1]);<br>
X=mexNormalize(X);<br>
Y=<span class="c002">randn</span>(100,1);<br>
Y=Y-repmat(<span class="c002">mean</span>(Y),[<span class="c002">size</span>(Y,1) 1]);<br>
Y=mexNormalize(Y);<br>
W0=<span class="c002">zeros</span>(<span class="c002">size</span>(X,2),<span class="c002">size</span>(Y,2));<br>
<span class="c001">% Regression experiments <br>
% 100 regression problems with the same design matrix X.</span><br>
<span class="c002">fprintf</span>(<span class="c003">'\nVarious regression experiments\n'</span>);<br>
param.compute_gram=true;<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Regression l1\n'</span>);<br>
param.loss=<span class="c003">'square'</span>;<br>
param.regul=<span class="c003">'l1'</span>;<br>
<span class="c001">% param.regul='group-lasso-l2';<br>
% param.size_group=10;</span><br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaFlat(Y,X,W0,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
param.regul=<span class="c003">'l1'</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'\nISTA + Regression l1\n'</span>);<br>
param.ista=true;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaFlat(Y,X,W0,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
param.regul=<span class="c003">'l1'</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'\nISTA + Regression l1 + Barzilai Borwein (SPARSA)\n'</span>);<br>
param.ista=true;<br>
param.barzilaiborwein=true;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaFlat(Y,X,W0,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
param.barzilaiborwein=false;<br>
<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\nSubgradient Descent + Regression l1\n'</span>);<br>
param.ista=false;<br>
param.subgrad=true;<br>
param.a=0.1;<br>
param.b=1000; <span class="c001">% arbitrary parameters</span><br>
max_it=param.max_it;<br>
it0=param.it0;<br>
param.max_it=500;<br>
param.it0=50;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaFlat(Y,X,W0,param);<br>
t=<span class="c002">toc</span>;<br>
param.subgrad=false;<br>
param.max_it=max_it;<br>
param.it0=it0;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Regression l2\n'</span>);<br>
param.regul=<span class="c003">'l2'</span>;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaFlat(Y,X,W0,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Regression l2 + sparse feature matrix\n'</span>);<br>
param.regul=<span class="c003">'l2'</span>;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaFlat(Y,<span class="c002">sparse</span>(X),W0,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Regression Elastic-Net\n'</span>);<br>
param.regul=<span class="c003">'elastic-net'</span>;<br>
param.lambda2=0.1;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaFlat(Y,X,W0,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Group Lasso L2\n'</span>);<br>
param.regul=<span class="c003">'group-lasso-l2'</span>;<br>
param.size_group=2;  <span class="c001">% all the groups are of size 2</span><br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaFlat(Y,X,W0,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Group Lasso L2 with variable size of groups \n'</span>);<br>
param.regul=<span class="c003">'group-lasso-l2'</span>;<br>
param2=param;<br>
param2.groups=int32(randi(5,1,<span class="c002">size</span>(X,2)));  <span class="c001">% all the groups are of size 2</span><br>
param2.lambda=10*param2.lambda;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaFlat(Y,X,W0,param2);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Trace Norm\n'</span>);<br>
param.regul=<span class="c003">'trace-norm-vec'</span>;<br>
param.size_group=5;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaFlat(Y,X,W0,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Regression Fused-Lasso\n'</span>);<br>
param.regul=<span class="c003">'fused-lasso'</span>;<br>
param.lambda2=0.1;<br>
param.lambda3=0.1; <span class="c001">%</span><br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaFlat(Y,X,W0,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Regression no regularization\n'</span>);<br>
param.regul=<span class="c003">'none'</span>;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaFlat(Y,X,W0,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Regression l1 with intercept \n'</span>);<br>
param.intercept=true;<br>
param.regul=<span class="c003">'l1'</span>;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaFlat(Y,[X ones(<span class="c002">size</span>(X,1),1)],[W0; <span class="c002">zeros</span>(1,<span class="c002">size</span>(W0,2))],param); <span class="c001">% adds a column of ones to X for the intercept</span><br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Regression l1 with intercept+ non-negative \n'</span>);<br>
param.pos=true;<br>
param.regul=<span class="c003">'l1'</span>;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaFlat(Y,[X ones(<span class="c002">size</span>(X,1),1)],[W0; <span class="c002">zeros</span>(1,<span class="c002">size</span>(W0,2))],param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
param.pos=false;<br>
param.intercept=false;<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\nISTA + Regression l0\n'</span>);<br>
param.regul=<span class="c003">'l0'</span>;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaFlat(Y,X,W0,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\nOne classification experiment\n'</span>);<br>
Y=2*double(<span class="c002">randn</span>(100,1) &gt; 0)-1;<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Logistic l1\n'</span>);<br>
param.regul=<span class="c003">'l1'</span>;<br>
param.loss=<span class="c003">'logistic'</span>;<br>
param.lambda=0.01;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaFlat(Y,X,W0,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<span class="c001">% can be used of course with other regularization functions, intercept,...</span><br>
param.regul=<span class="c003">'l1'</span>;<br>
param.loss=<span class="c003">'weighted-logistic'</span>;<br>
param.lambda=0.01;<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + weighted Logistic l1 + sparse matrix\n'</span>);<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaFlat(Y,X,W0,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c001">% can be used of course with other regularization functions, intercept,...</span><br>
<br>
param.loss=<span class="c003">'logistic'</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Logistic l1 + sparse matrix\n'</span>);<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaFlat(Y,<span class="c002">sparse</span>(X),W0,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<span class="c001">% can be used of course with other regularization functions, intercept,...<br>
<br>
% Multi-Class classification</span><br>
Y=double(<span class="c002">ceil</span>(5*<span class="c002">rand</span>(100,1000))-1);<br>
param.loss=<span class="c003">'multi-logistic'</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Multi-Class Logistic l1\n'</span>);<br>
<span class="c002">tic</span><br>
nclasses=<span class="c002">max</span>(Y(:))+1;<br>
W0=<span class="c002">zeros</span>(<span class="c002">size</span>(X,2),nclasses*<span class="c002">size</span>(Y,2));<br>
[W optim_info]=mexFistaFlat(Y,X,W0,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<span class="c002">pause</span><br>
<span class="c001">% can be used of course with other regularization functions, intercept,...<br>
<br>
% Multi-Task regression</span><br>
Y=<span class="c002">randn</span>(100,100);<br>
Y=Y-repmat(<span class="c002">mean</span>(Y),[<span class="c002">size</span>(Y,1) 1]);<br>
Y=mexNormalize(Y);<br>
param.compute_gram=false;<br>
W0=<span class="c002">zeros</span>(<span class="c002">size</span>(X,2),<span class="c002">size</span>(Y,2));<br>
param.loss=<span class="c003">'square'</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Regression l1l2 \n'</span>);<br>
param.regul=<span class="c003">'l1l2'</span>;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaFlat(Y,X,W0,param);<br>
<span class="c002">toc<br>
fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Regression l1linf \n'</span>);<br>
param.regul=<span class="c003">'l1linf'</span>;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaFlat(Y,X,W0,param);<br>
<span class="c002">toc<br>
fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Regression l1l2 + l1 \n'</span>);<br>
param.regul=<span class="c003">'l1l2+l1'</span>;<br>
param.lambda2=0.1;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaFlat(Y,X,W0,param);<br>
<span class="c002">toc<br>
fprintf</span>(<span class="c003">'mean loss: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Regression l1linf + l1 \n'</span>);<br>
param.regul=<span class="c003">'l1linf+l1'</span>;<br>
param.lambda2=0.1;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaFlat(Y,X,W0,param);<br>
<span class="c002">toc<br>
fprintf</span>(<span class="c003">'mean loss: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Regression l1linf + row + columns \n'</span>);<br>
param.regul=<span class="c003">'l1linf-row-column'</span>;<br>
param.lambda2=0.1;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaFlat(Y,X,W0,param);<br>
<span class="c002">toc<br>
fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
<span class="c001">% Multi-Task Classification</span><br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Logistic + l1l2 \n'</span>);<br>
param.regul=<span class="c003">'l1l2'</span>;<br>
param.loss=<span class="c003">'logistic'</span>;<br>
Y=2*double(<span class="c002">randn</span>(100,100) &gt; 0)-1;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaFlat(Y,X,W0,param);<br>
<span class="c002">toc<br>
fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<span class="c001">% Multi-Class + Multi-Task Regularization</span><br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Multi-Class Logistic l1l2 \n'</span>);<br>
Y=double(<span class="c002">ceil</span>(5*<span class="c002">rand</span>(100,1000))-1);<br>
param.loss=<span class="c003">'multi-logistic'</span>;<br>
param.regul=<span class="c003">'l1l2'</span>;<br>
<span class="c002">tic</span><br>
nclasses=<span class="c002">max</span>(Y(:))+1;<br>
W0=<span class="c002">zeros</span>(<span class="c002">size</span>(X,2),nclasses*<span class="c002">size</span>(Y,2));<br>
[W optim_info]=mexFistaFlat(Y,X,W0,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<span class="c001">% can be used of course with other regularization functions, intercept,...</span></td></tr>
</table>
<h3 class="subsection" id="sec38">5.9  Function mexFistaTree</h3>
<p>
Given a matrix <span class="c009">X</span>=[<span class="c009">x</span><sup>1</sup>,…,<span class="c009">x</span><sup><span class="c007">p</span></sup>]<sup><span class="c007">T</span></sup> in ℝ<sup><span class="c007">m</span> × <span class="c007">p</span></sup>, and a matrix <span class="c009">Y</span>=[<span class="c009">y</span><sup>1</sup>,…,<span class="c009">y</span><sup><span class="c007">n</span></sup>], it solves the optimization problems presented in the previous section, with the same regularization functions as mexProximalTree.
see usage details below:</p><table class="lstframe c011"><tr><td class="mouselstlisting"><span class="c001">% <br>
% Usage: [W [optim]]=mexFistaTree(Y,X,W0,tree,param);<br>
%<br>
% Name: mexFistaTree<br>
%<br>
% Description: mexFistaTree solves sparse regularized problems.<br>
%         X is a design matrix of size m x p<br>
%         X=[x^1,...,x^n]', where the x_i's are the rows of X<br>
%         Y=[y^1,...,y^n] is a matrix of size m x n<br>
%         It implements the algorithms FISTA, ISTA and subgradient descent for solving<br>
%<br>
%           min_W  loss(W) + lambda psi(W)<br>
%          <br>
%         The function psi are those used by mexProximalTree (see documentation)<br>
%         for the loss functions, see the documentation of mexFistaFlat<br>
%<br>
%         This function can also handle intercepts (last row of W is not regularized),<br>
%         and/or non-negativity constraints on W and sparse matrices X<br>
%<br>
% Inputs: Y:  double dense m x n matrix<br>
%         X:  double dense or sparse m x p matrix   <br>
%         W0:  double dense p x n matrix or p x Nn matrix (for multi-logistic loss)<br>
%              initial guess<br>
%         tree: struct (see documentation of mexProximalTree)<br>
%         param: struct<br>
%            param.loss (choice of loss, see above)<br>
%            param.regul (choice of regularization, see function mexProximalFlat)<br>
%            param.lambda (regularization parameter)<br>
%            param.lambda2 (optional, regularization parameter, 0 by default)<br>
%            param.lambda3 (optional, regularization parameter, 0 by default)<br>
%            param.verbose (optional, verbosity level, false by default)<br>
%            param.pos (optional, adds positivity constraints on the<br>
%                coefficients, false by default)<br>
%            param.transpose (optional, transpose the matrix in the regularization function)<br>
%            param.size_group (optional, for regularization functions assuming a group<br>
%                 structure)<br>
%            param.numThreads (optional, number of threads for exploiting<br>
%                multi-core / multi-cpus. By default, it takes the value -1,<br>
%                which automatically selects all the available CPUs/cores).<br>
%            param.max_it (optional, maximum number of iterations, 100 by default)<br>
%            param.it0 (optional, frequency for computing duality gap, every 10 iterations by default)<br>
%            param.tol (optional, tolerance for stopping criteration, which is a relative duality gap<br>
%                if it is available, or a relative change of parameters).<br>
%            param.gamma (optional, multiplier for increasing the parameter L in fista, 1.5 by default)<br>
%            param.L0 (optional, initial parameter L in fista, 0.1 by default, should be small enough)<br>
%            param.fixed_step (deactive the line search for L in fista and use param.L0 instead)<br>
%            param.compute_gram (optional, pre-compute X^TX, false by default).<br>
%            param.intercept (optional, do not regularize last row of W, false by default).<br>
%            param.ista (optional, use ista instead of fista, false by default).<br>
%            param.subgrad (optional, if not param.ista, use subradient descent instead of fista, false by default).<br>
%            param.a, param.b (optional, if param.subgrad, the gradient step is a/(t+b)<br>
%            also similar options as mexProximalTree<br>
%<br>
%            the function also implements the ADMM algorithm via an option param.admm=true. It is not documented<br>
%            and you need to look at the source code to use it.<br>
%<br>
% Output:  W:  double dense p x n matrix or p x Nn matrix (for multi-logistic loss)<br>
%          optim: optional, double dense 4 x n matrix.<br>
%              first row: values of the objective functions.<br>
%              third row: values of the relative duality gap (if available)<br>
%              fourth row: number of iterations<br>
%<br>
% Author: Julien Mairal, 2010</span></td></tr>
</table><p>
The following piece of code illustrates how to use this function.
</p><table class="lstframe c011"><tr><td class="mouselstlisting"><span class="c002">format</span> compact;<br>
param.num_threads=-1; <span class="c001">% all cores (-1 by default)</span><br>
param.verbose=false;   <span class="c001">% verbosity, false by default</span><br>
param.lambda=0.001; <span class="c001">% regularization parameter</span><br>
param.it0=10;      <span class="c001">% frequency for duality gap computations</span><br>
param.max_it=200; <span class="c001">% maximum number of iterations</span><br>
param.L0=0.1;<br>
param.tol=1e-5;<br>
param.intercept=false;<br>
param.pos=false;<br>
<br>
<span class="c001">% Example 2 of tree structure<br>
% tree structured groups:<br>
% g1= {0 1 2 3 4 5 6 7 8 9}    root(g1) = { };<br>
% g2= {0 1 2 3 4 5}            root(g2) = {0 1 2};<br>
% g3= {3 4}                    root(g3) = {3 4};<br>
% g4= {5}                      root(g4) = {5};<br>
% g5= {6 7 8 9}                root(g5) = { };<br>
% g6= {6 7}                    root(g6) = {6 7};<br>
% g7= {8 9}                    root(g7) = {8};<br>
% g8 = {9}                     root(g8) = {9};</span><br>
tree.own_variables=  int32([0 0 3 5 6 6 8 9]);   <span class="c001">% pointer to the first variable of each group</span><br>
tree.N_own_variables=int32([0 3 2 1 0 2 1 1]); <span class="c001">% number of "root" variables in each group</span><br>
tree.eta_g=[1 1 1 2 2 2 2.5 2.5];<br>
tree.groups=<span class="c002">sparse</span>([0 0 0 0 0 0 0 0; ...<br>
                    1 0 0 0 0 0 0 0; ...<br>
                    0 1 0 0 0 0 0 0; ...<br>
                    0 1 0 0 0 0 0 0; ...<br>
                    1 0 0 0 0 0 0 0; ...<br>
                    0 0 0 0 1 0 0 0; ...<br>
                    0 0 0 0 1 0 0 0; ...<br>
                    0 0 0 0 0 0 1 0]);  <span class="c001">% first group should always be the root of the tree</span><br>
<br>
X=<span class="c002">randn</span>(100,10);<br>
X=X-repmat(<span class="c002">mean</span>(X),[<span class="c002">size</span>(X,1) 1]);<br>
X=mexNormalize(X);<br>
Y=<span class="c002">randn</span>(100,100);<br>
Y=Y-repmat(<span class="c002">mean</span>(Y),[<span class="c002">size</span>(Y,1) 1]);<br>
Y=mexNormalize(Y);<br>
W0=<span class="c002">zeros</span>(<span class="c002">size</span>(X,2),<span class="c002">size</span>(Y,2));<br>
<span class="c001">% Regression experiments <br>
% 100 regression problems with the same design matrix X.</span><br>
<span class="c002">fprintf</span>(<span class="c003">'\nVarious regression experiments\n'</span>);<br>
param.compute_gram=true;<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Regression tree-l2\n'</span>);<br>
param.loss=<span class="c003">'square'</span>;<br>
param.regul=<span class="c003">'tree-l2'</span>;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaTree(Y,X,W0,tree,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Regression tree-linf\n'</span>);<br>
param.regul=<span class="c003">'tree-linf'</span>;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaTree(Y,X,W0,tree,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
<span class="c001">% works also with non tree-structured regularization. tree is ignored</span><br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Regression Fused-Lasso\n'</span>);<br>
param.regul=<span class="c003">'fused-lasso'</span>;<br>
param.lambda2=0.001;<br>
param.lambda3=0.001; <span class="c001">%</span><br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaTree(Y,X,W0,tree,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\nISTA + Regression tree-l0\n'</span>);<br>
param.regul=<span class="c003">'tree-l0'</span>;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaTree(Y,X,W0,tree,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Regression tree-l2 with intercept \n'</span>);<br>
param.intercept=true;<br>
param.regul=<span class="c003">'tree-l2'</span>;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaTree(Y,[X ones(<span class="c002">size</span>(X,1),1)],[W0; <span class="c002">zeros</span>(1,<span class="c002">size</span>(W0,2))],tree,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
param.intercept=false;<br>
<br>
<span class="c001">% Classification</span><br>
<span class="c002">fprintf</span>(<span class="c003">'\nOne classification experiment\n'</span>);<br>
Y=2*double(<span class="c002">randn</span>(100,<span class="c002">size</span>(Y,2)) &gt; 0)-1;<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Logistic + tree-linf\n'</span>);<br>
param.regul=<span class="c003">'tree-linf'</span>;<br>
param.loss=<span class="c003">'logistic'</span>;<br>
param.lambda=0.001;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaTree(Y,X,W0,tree,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<span class="c001">% can be used of course with other regularization functions, intercept,...<br>
<br>
% Multi-Class classification</span><br>
Y=double(<span class="c002">ceil</span>(5*<span class="c002">rand</span>(100,<span class="c002">size</span>(Y,2)))-1);<br>
param.loss=<span class="c003">'multi-logistic'</span>;<br>
param.regul=<span class="c003">'tree-l2'</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Multi-Class Logistic + tree-l2 \n'</span>);<br>
<span class="c002">tic</span><br>
nclasses=<span class="c002">max</span>(Y(:))+1;<br>
W0=<span class="c002">zeros</span>(<span class="c002">size</span>(X,2),nclasses*<span class="c002">size</span>(Y,2));<br>
[W optim_info]=mexFistaTree(Y,X,W0,tree,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<span class="c001">% can be used of course with other regularization functions, intercept,...<br>
<br>
% Multi-Task regression</span><br>
Y=<span class="c002">randn</span>(100,<span class="c002">size</span>(Y,2));<br>
Y=Y-repmat(<span class="c002">mean</span>(Y),[<span class="c002">size</span>(Y,1) 1]);<br>
Y=mexNormalize(Y);<br>
param.compute_gram=false;<br>
param.verbose=true;   <span class="c001">% verbosity, false by default</span><br>
W0=<span class="c002">zeros</span>(<span class="c002">size</span>(X,2),<span class="c002">size</span>(Y,2));<br>
param.loss=<span class="c003">'square'</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Regression multi-task-tree \n'</span>);<br>
param.regul=<span class="c003">'multi-task-tree'</span>;<br>
param.lambda2=0.001;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaTree(Y,X,W0,tree,param);<br>
<span class="c002">toc<br>
fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
<span class="c001">% Multi-Task Classification</span><br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Logistic + multi-task-tree \n'</span>);<br>
param.regul=<span class="c003">'multi-task-tree'</span>;<br>
param.lambda2=0.001;<br>
param.loss=<span class="c003">'logistic'</span>;<br>
Y=2*double(<span class="c002">randn</span>(100,<span class="c002">size</span>(Y,2)) &gt; 0)-1;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaTree(Y,X,W0,tree,param);<br>
<span class="c002">toc<br>
fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
<span class="c001">% Multi-Class + Multi-Task Regularization</span><br>
param.verbose=false;<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Multi-Class Logistic +multi-task-tree \n'</span>);<br>
Y=double(<span class="c002">ceil</span>(5*<span class="c002">rand</span>(100,<span class="c002">size</span>(Y,2)))-1);<br>
param.loss=<span class="c003">'multi-logistic'</span>;<br>
param.regul=<span class="c003">'multi-task-tree'</span>;<br>
<span class="c002">tic</span><br>
nclasses=<span class="c002">max</span>(Y(:))+1;<br>
W0=<span class="c002">zeros</span>(<span class="c002">size</span>(X,2),nclasses*<span class="c002">size</span>(Y,2));<br>
[W optim_info]=mexFistaTree(Y,X,W0,tree,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<span class="c001">% can be used of course with other regularization functions, intercept,...</span><br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Multi-Class Logistic +multi-task-tree + sparse matrix \n'</span>);<br>
<span class="c002">tic</span><br>
nclasses=<span class="c002">max</span>(Y(:))+1;<br>
W0=<span class="c002">zeros</span>(<span class="c002">size</span>(X,2),nclasses*<span class="c002">size</span>(Y,2));<br>
[W optim_info]=mexFistaTree(Y,<span class="c002">sparse</span>(X),W0,tree,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));</td></tr>
</table>
<h3 class="subsection" id="sec39">5.10  Function mexFistaGraph</h3>
<p>
Given a matrix <span class="c009">X</span>=[<span class="c009">x</span><sup>1</sup>,…,<span class="c009">x</span><sup><span class="c007">p</span></sup>]<sup><span class="c007">T</span></sup> in ℝ<sup><span class="c007">m</span> × <span class="c007">p</span></sup>, and a matrix <span class="c009">Y</span>=[<span class="c009">y</span><sup>1</sup>,…,<span class="c009">y</span><sup><span class="c007">n</span></sup>], it solves the optimization problems presented in the previous section, with the same regularization functions as mexProximalGraph.
see usage details below:</p><table class="lstframe c011"><tr><td class="mouselstlisting"><span class="c001">% <br>
% Usage: [W [optim]]=mexFistaGraph(Y,X,W0,graph,param);<br>
%<br>
% Name: mexFistaGraph<br>
%<br>
% Description: mexFistaGraph solves sparse regularized problems.<br>
%         X is a design matrix of size m x p<br>
%         X=[x^1,...,x^n]', where the x_i's are the rows of X<br>
%         Y=[y^1,...,y^n] is a matrix of size m x n<br>
%         It implements the algorithms FISTA, ISTA and subgradient descent.<br>
%<br>
%         It implements the algorithms FISTA, ISTA and subgradient descent for solving<br>
%<br>
%           min_W  loss(W) + lambda psi(W)<br>
%          <br>
%         The function psi are those used by mexProximalGraph (see documentation)<br>
%         for the loss functions, see the documentation of mexFistaFlat<br>
%         <br>
%         This function can also handle intercepts (last row of W is not regularized),<br>
%         and/or non-negativity constraints on W.<br>
%<br>
% Inputs: Y:  double dense m x n matrix<br>
%         X:  double dense or sparse m x p matrix   <br>
%         W0:  double dense p x n matrix or p x Nn matrix (for multi-logistic loss)<br>
%              initial guess<br>
%         graph: struct (see documentation of mexProximalGraph)<br>
%         param: struct<br>
%            param.loss (choice of loss, see above)<br>
%            param.regul (choice of regularization, see function mexProximalFlat)<br>
%            param.lambda (regularization parameter)<br>
%            param.lambda2 (optional, regularization parameter, 0 by default)<br>
%            param.lambda3 (optional, regularization parameter, 0 by default)<br>
%            param.verbose (optional, verbosity level, false by default)<br>
%            param.pos (optional, adds positivity constraints on the<br>
%                coefficients, false by default)<br>
%            param.numThreads (optional, number of threads for exploiting<br>
%                multi-core / multi-cpus. By default, it takes the value -1,<br>
%                which automatically selects all the available CPUs/cores).<br>
%            param.max_it (optional, maximum number of iterations, 100 by default)<br>
%            param.it0 (optional, frequency for computing duality gap, every 10 iterations by default)<br>
%            param.tol (optional, tolerance for stopping criteration, which is a relative duality gap<br>
%                if it is available, or a relative change of parameters).<br>
%            param.gamma (optional, multiplier for increasing the parameter L in fista, 1.5 by default)<br>
%            param.L0 (optional, initial parameter L in fista, 0.1 by default, should be small enough)<br>
%            param.fixed_step (deactive the line search for L in fista and use param.L0 instead)<br>
%            param.compute_gram (optional, pre-compute X^TX, false by default).<br>
%            param.intercept (optional, do not regularize last row of W, false by default).<br>
%            param.ista (optional, use ista instead of fista, false by default).<br>
%            param.subgrad (optional, if not param.ista, use subradient descent instead of fista, false by default).<br>
%            param.a, param.b (optional, if param.subgrad, the gradient step is a/(t+b)<br>
%            also similar options as mexProximalTree<br>
%<br>
%            the function also implements the ADMM algorithm via an option param.admm=true. It is not documented<br>
%            and you need to look at the source code to use it.<br>
%<br>
%<br>
% Output:  W:  double dense p x n matrix or p x Nn matrix (for multi-logistic loss)<br>
%          optim: optional, double dense 4 x n matrix.<br>
%              first row: values of the objective functions.<br>
%              third row: values of the relative duality gap (if available)<br>
%              fourth row: number of iterations<br>
%<br>
% Author: Julien Mairal, 2010</span></td></tr>
</table><p>
The following piece of code illustrates how to use this function.
</p><table class="lstframe c011"><tr><td class="mouselstlisting"><span class="c002">clear</span> <span class="c002">all</span>;<br>
<span class="c002">randn</span>(<span class="c003">'seed'</span>,0);<br>
<span class="c002">format</span> compact;<br>
param.num_threads=-1; <span class="c001">% all cores (-1 by default)</span><br>
param.verbose=false;   <span class="c001">% verbosity, false by default</span><br>
param.lambda=0.1; <span class="c001">% regularization parameter</span><br>
param.it0=1;      <span class="c001">% frequency for duality gap computations</span><br>
param.max_it=100; <span class="c001">% maximum number of iterations</span><br>
param.L0=0.1;<br>
param.tol=1e-5;<br>
param.intercept=false;<br>
param.pos=false;<br>
<br>
graph.eta_g=[1 1 1 1 1];<br>
graph.groups=<span class="c002">sparse</span>([0 0 0 1 0;<br>
                     0 0 0 0 0;<br>
                     0 0 0 0 0;<br>
                     0 0 0 0 0;<br>
                     0 0 1 0 0]);   <span class="c001">% g5 is included in g3, and g2 is included in g4</span><br>
graph.groups_var=<span class="c002">sparse</span>([1 0 0 0 0;<br>
                         1 0 0 0 0;<br>
                         1 0 0 0 0 ;<br>
                         1 1 0 0 0;<br>
                         0 1 0 1 0;<br>
                         0 1 0 1 0;<br>
                         0 1 0 0 1;<br>
                         0 0 0 0 1;<br>
                         0 0 0 0 1;<br>
                         0 0 1 0 0]); <span class="c001">% represents direct inclusion relations </span><br>
<br>
X=<span class="c002">randn</span>(100,10);<br>
param.verbose=true;<br>
<span class="c001">%X=eye(10);</span><br>
X=X-repmat(<span class="c002">mean</span>(X),[<span class="c002">size</span>(X,1) 1]);<br>
X=mexNormalize(X);<br>
Y=<span class="c002">randn</span>(100,1);<br>
Y=Y-repmat(<span class="c002">mean</span>(Y),[<span class="c002">size</span>(Y,1) 1]);<br>
Y=mexNormalize(Y);<br>
W0=<span class="c002">zeros</span>(<span class="c002">size</span>(X,2),<span class="c002">size</span>(Y,2));<br>
<span class="c001">% Regression experiments <br>
% 100 regression problems with the same design matrix X.</span><br>
<span class="c002">fprintf</span>(<span class="c003">'\nVarious regression experiments\n'</span>);<br>
param.compute_gram=true;<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Regression graph\n'</span>);<br>
param.loss=<span class="c003">'square'</span>;<br>
param.regul=<span class="c003">'graph'</span>;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaGraph(Y,X,W0,graph,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\nADMM + Regression graph\n'</span>);<br>
param.admm=true;<br>
param.lin_admm=true;<br>
param.c=1;<br>
param.delta=1;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaGraph(Y,X,W0,graph,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, stopping criterion: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
param.admm=false;<br>
param.max_it=5;<br>
param.it0=1;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaGraph(Y,X,W,graph,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
<span class="c001">% works also with non graph-structured regularization. graph is ignored</span><br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Regression Fused-Lasso\n'</span>);<br>
param.regul=<span class="c003">'fused-lasso'</span>;<br>
param.lambda2=0.01;<br>
param.lambda3=0.01; <span class="c001">%</span><br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaGraph(Y,X,W0,graph,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Regression graph with intercept \n'</span>);<br>
param.intercept=true;<br>
param.regul=<span class="c003">'graph'</span>;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaGraph(Y,[X ones(<span class="c002">size</span>(X,1),1)],[W0; <span class="c002">zeros</span>(1,<span class="c002">size</span>(W0,2))],graph,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
param.intercept=false;<br>
<br>
<span class="c001">% Classification</span><br>
<span class="c002">fprintf</span>(<span class="c003">'\nOne classification experiment\n'</span>);<br>
Y=2*double(<span class="c002">randn</span>(100,<span class="c002">size</span>(Y,2)) &gt; 0)-1;<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Logistic + graph-linf\n'</span>);<br>
param.regul=<span class="c003">'graph'</span>;<br>
param.loss=<span class="c003">'logistic'</span>;<br>
param.lambda=0.01;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaGraph(Y,X,W0,graph,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<span class="c001">% can be used of course with other regularization functions, intercept,...<br>
<br>
% Multi-Class classification</span><br>
Y=double(<span class="c002">ceil</span>(5*<span class="c002">rand</span>(100,<span class="c002">size</span>(Y,2)))-1);<br>
param.loss=<span class="c003">'multi-logistic'</span>;<br>
param.regul=<span class="c003">'graph'</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Multi-Class Logistic + graph \n'</span>);<br>
<span class="c002">tic</span><br>
nclasses=<span class="c002">max</span>(Y(:))+1;<br>
W0=<span class="c002">zeros</span>(<span class="c002">size</span>(X,2),nclasses*<span class="c002">size</span>(Y,2));<br>
[W optim_info]=mexFistaGraph(Y,X,W0,graph,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<span class="c001">% can be used of course with other regularization functions, intercept,...<br>
<br>
% Multi-Task regression</span><br>
Y=<span class="c002">randn</span>(100,<span class="c002">size</span>(Y,2));<br>
Y=Y-repmat(<span class="c002">mean</span>(Y),[<span class="c002">size</span>(Y,1) 1]);<br>
Y=mexNormalize(Y);<br>
param.compute_gram=false;<br>
param.verbose=true;   <span class="c001">% verbosity, false by default</span><br>
W0=<span class="c002">zeros</span>(<span class="c002">size</span>(X,2),<span class="c002">size</span>(Y,2));<br>
param.loss=<span class="c003">'square'</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Regression multi-task-graph \n'</span>);<br>
param.regul=<span class="c003">'multi-task-graph'</span>;<br>
param.lambda2=0.01;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaGraph(Y,X,W0,graph,param);<br>
<span class="c002">toc<br>
fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<br>
<span class="c001">% Multi-Task Classification</span><br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Logistic + multi-task-graph \n'</span>);<br>
param.regul=<span class="c003">'multi-task-graph'</span>;<br>
param.lambda2=0.01;<br>
param.loss=<span class="c003">'logistic'</span>;<br>
Y=2*double(<span class="c002">randn</span>(100,<span class="c002">size</span>(Y,2)) &gt; 0)-1;<br>
<span class="c002">tic</span><br>
[W optim_info]=mexFistaGraph(Y,X,W0,graph,param);<br>
<span class="c002">toc<br>
fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<span class="c001">% Multi-Class + Multi-Task Regularization</span><br>
<br>
param.verbose=false;<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Multi-Class Logistic +multi-task-graph \n'</span>);<br>
Y=double(<span class="c002">ceil</span>(5*<span class="c002">rand</span>(100,<span class="c002">size</span>(Y,2)))-1);<br>
param.loss=<span class="c003">'multi-logistic'</span>;<br>
param.regul=<span class="c003">'multi-task-graph'</span>;<br>
<span class="c002">tic</span><br>
nclasses=<span class="c002">max</span>(Y(:))+1;<br>
W0=<span class="c002">zeros</span>(<span class="c002">size</span>(X,2),nclasses*<span class="c002">size</span>(Y,2));<br>
[W optim_info]=mexFistaGraph(Y,X,W0,graph,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
<span class="c001">% can be used of course with other regularization functions, intercept,...</span></td></tr>
</table>
<h3 class="subsection" id="sec40">5.11  Function mexFistaPathCoding</h3>
<p>
Similarly, the toolbox handles the penalties of [<a href="doc_spams010.html#mairal14">24</a>] with the following function
</p><table class="lstframe c011"><tr><td class="mouselstlisting"><span class="c001">% <br>
% Usage: [W [optim]]=mexFistaPathCoding(Y,X,W0,DAG,param);<br>
%<br>
% Name: mexFistaPathCoding<br>
%<br>
% Description: mexFistaPathCoding solves sparse regularized problems for the <br>
%         path coding penalties of http://arxiv.org/abs/1204.4539<br>
%         X is a design matrix of size m x p<br>
%         X=[x^1,...,x^n]', where the x_i's are the rows of X<br>
%         Y=[y^1,...,y^n] is a matrix of size m x n<br>
%         It implements the algorithms FISTA, ISTA and subgradient descent.<br>
%<br>
%         It implements the algorithms FISTA, ISTA and subgradient descent for solving<br>
%<br>
%           min_W  loss(W) + lambda psi(W)<br>
%          <br>
%         The function psi are those used by mexProximalPathCoding (see documentation)<br>
%         for the loss functions, see the documentation of mexFistaFlat<br>
%         <br>
%         This function can also handle intercepts (last row of W is not regularized),<br>
%         and/or non-negativity constraints on W.<br>
%<br>
% Inputs: Y:  double dense m x n matrix<br>
%         X:  double dense or sparse m x p matrix   <br>
%         W0:  double dense p x n matrix or p x Nn matrix (for multi-logistic loss)<br>
%              initial guess<br>
%         DAG: struct (see documentation of mexProximalPathCoding)<br>
%         param: struct<br>
%            param.loss (choice of loss, see above)<br>
%            param.regul (choice of regularization, see function mexProximalPathCoding)<br>
%            param.lambda (regularization parameter)<br>
%            param.lambda2 (optional, regularization parameter, 0 by default)<br>
%            param.lambda3 (optional, regularization parameter, 0 by default)<br>
%            param.verbose (optional, verbosity level, false by default)<br>
%            param.pos (optional, adds positivity constraints on the<br>
%                coefficients, false by default)<br>
%            param.numThreads (optional, number of threads for exploiting<br>
%                multi-core / multi-cpus. By default, it takes the value -1,<br>
%                which automatically selects all the available CPUs/cores).<br>
%            param.max_it (optional, maximum number of iterations, 100 by default)<br>
%            param.it0 (optional, frequency for computing duality gap, every 10 iterations by default)<br>
%            param.tol (optional, tolerance for stopping criteration, which is a relative duality gap<br>
%                if it is available, or a relative change of parameters).<br>
%            param.gamma (optional, multiplier for increasing the parameter L in fista, 1.5 by default)<br>
%            param.L0 (optional, initial parameter L in fista, 0.1 by default, should be small enough)<br>
%            param.fixed_step (deactive the line search for L in fista and use param.L0 instead)<br>
%            param.compute_gram (optional, pre-compute X^TX, false by default).<br>
%            param.intercept (optional, do not regularize last row of W, false by default).<br>
%            param.ista (optional, use ista instead of fista, false by default).<br>
%            param.subgrad (optional, if not param.ista, use subradient descent instead of fista, false by default).<br>
%            param.a, param.b (optional, if param.subgrad, the gradient step is a/(t+b)<br>
%            also similar options as mexProximalPathCoding<br>
%<br>
%<br>
% Output:  W:  double dense p x n matrix or p x Nn matrix (for multi-logistic loss)<br>
%          optim: optional, double dense 4 x n matrix.<br>
%              first row: values of the objective functions.<br>
%              third row: values of the relative duality gap (if available)<br>
%              fourth row: number of iterations<br>
%<br>
% Author: Julien Mairal, 2012</span></td></tr>
</table><p>
The following piece of code illustrates how to use this function.
</p><table class="lstframe c011"><tr><td class="mouselstlisting"><span class="c002">clear</span> <span class="c002">all</span>;<br>
<span class="c002">rand</span>(<span class="c003">'seed'</span>,0);<br>
<span class="c002">randn</span>(<span class="c003">'seed'</span>,0);<br>
<span class="c002">fprintf</span>(<span class="c003">'test mexFistaPathCoding\n'</span>);<br>
p=100;<br>
n=1000;<br>
<span class="c001">% generate a DAG</span><br>
G=sprand(p,p,0.05);<br>
G=mexRemoveCyclesGraph(G);<br>
<span class="c002">fprintf</span>(<span class="c003">'\n'</span>);<br>
<br>
<span class="c001">% generate a data matrix</span><br>
X=<span class="c002">randn</span>(n,p);<br>
X=X-repmat(<span class="c002">mean</span>(X),[<span class="c002">size</span>(X,1) 1]);<br>
X=mexNormalize(X);<br>
Y=<span class="c002">randn</span>(n,2);<br>
Y=Y-repmat(<span class="c002">mean</span>(Y),[<span class="c002">size</span>(Y,1) 1]);<br>
Y=mexNormalize(Y);<br>
W0=<span class="c002">zeros</span>(<span class="c002">size</span>(X,2),<span class="c002">size</span>(Y,2));<br>
<br>
<span class="c001">% input graph</span><br>
graph.weights=G;<br>
graph.stop_weights=<span class="c002">zeros</span>(1,p);<br>
graph.start_weights=10*ones(1,p);<br>
<br>
<span class="c001">% FISTA parameters</span><br>
param.num_threads=-1; <span class="c001">% all cores (-1 by default)</span><br>
param.verbose=true;   <span class="c001">% verbosity, false by default</span><br>
param.lambda=0.005; <span class="c001">% regularization parameter</span><br>
param.it0=1;      <span class="c001">% frequency for duality gap computations</span><br>
param.max_it=100; <span class="c001">% maximum number of iterations</span><br>
param.L0=0.01;<br>
param.tol=1e-4;<br>
param.precision=10000000;<br>
param.pos=false;<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'Square Loss + convex path penalty\n'</span>);<br>
param.loss=<span class="c003">'square'</span>;<br>
param.regul=<span class="c003">'graph-path-conv'</span>;<br>
<span class="c002">tic</span><br>
[W1 optim_info]=mexFistaPathCoding(Y,X,W0,graph,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));<br>
num=mexCountConnexComponents(graph.weights,W1(:,1));<br>
<span class="c002">fprintf</span>(<span class="c003">'Num of connected components: %d\n'</span>,num);<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\n'</span>);<br>
<span class="c002">fprintf</span>(<span class="c003">'Square Loss + non-convex path penalty\n'</span>);<br>
param.loss=<span class="c003">'square'</span>;<br>
param.regul=<span class="c003">'graph-path-l0'</span>;<br>
param.lambda=0.0001; <span class="c001">% regularization parameter</span><br>
param.ista=true;<br>
<span class="c002">tic</span><br>
[W2 optim_info]=mexFistaPathCoding(Y,X,W0,graph,param);<br>
t=<span class="c002">toc</span>;<br>
num=mexCountConnexComponents(graph.weights,W2(:,1));<br>
<span class="c002">fprintf</span>(<span class="c003">'Num of connected components: %d\n'</span>,num);<br>
<br>
<span class="c002">fprintf</span>(<span class="c003">'\n'</span>);<br>
<span class="c002">fprintf</span>(<span class="c003">'Note that for non-convex penalties, continuation strategies sometimes perform better:\n'</span>);<br>
tablambda=param.lambda*<span class="c002">sqrt</span>(<span class="c002">sqrt</span>(<span class="c002">sqrt</span>(2))).^(20:-1:0);<br>
lambda_orig=param.lambda;<br>
<span class="c002">tic</span><br>
W2=W0;<br>
<span class="c002">for</span> ii = 1:<span class="c002">length</span>(tablambda)<br>
   param.lambda=tablambda(ii);<br>
   param.verbose=false;<br>
   [W2]=mexFistaPathCoding(Y,X,W2,graph,param);<br>
<span class="c002">end</span><br>
param.verbose=true;<br>
param.lambda=lambda_orig;<br>
[W2 optim_info]=mexFistaPathCoding(Y,X,W2,graph,param);<br>
t=<span class="c002">toc</span>;<br>
num=mexCountConnexComponents(graph.weights,W2(:,1));<br>
param.ista=false;<br>
<span class="c002">fprintf</span>(<span class="c003">'Num of connected components: %d\n'</span>,num);</td></tr>
</table>
<h3 class="subsection" id="sec41">5.12  Function solverPoisson</h3>
<p>
The following problem is addressed here
</p><table class="display dcenter"><tr class="c016"><td class="dcell">   </td><td class="dcell"><table class="display"><tr><td class="dcell c012">&nbsp;</td></tr>
<tr><td class="dcell c012">min</td></tr>
<tr><td class="dcell c012"><span class="c009">w</span> ∈ ℝ<sub>+</sub><sup><span class="c007">p</span></sup></td></tr>
</table></td><td class="dcell">  </td><td class="dcell"><table class="display"><tr><td class="dcell c012"><span class="c007">p</span></td></tr>
<tr><td class="dcell c012"><span class="c006">∑</span></td></tr>
<tr><td class="dcell c012"><span class="c007">j</span>=1</td></tr>
</table></td><td class="dcell"> <span class="c009">x</span><sub><span class="c007">i</span></sub><sup>⊤</sup><span class="c009">w</span> + δ − <span class="c007">y</span><sub><span class="c007">i</span></sub> log(<span class="c009">x</span><sub><span class="c007">i</span></sub><sup>⊤</sup><span class="c009">w</span> +δ) + λ ψ(<span class="c009">w</span>).
</td></tr>
</table><table class="lstframe c011"><tr><td class="mouselstlisting"><span class="c001">% <br>
% Usage:   [W, [optim]]=solverPoisson(Y,X,W0,param);<br>
%<br>
% Name: solverPoisson<br>
%<br>
% Description: solverPoisson solves the regularized Poisson regression<br>
% problem for every column of Y:<br>
%<br>
%   min_w  \sum_i (x_i' w) + delta  - y_i log( x_i' w + delta) + lambda psi(w) (1)<br>
%<br>
% delta should be positive. The method solves (1) for decreasing values of delta,<br>
% until it reaches the desired target value. The algorithm ISTA is used<br>
% with Barzilai-Borwein steps.<br>
%<br>
% Inputs: Y:  double m x n matrix (non-negative values)<br>
%         X:  double m x p matrix (non-negative values)<br>
%         W0:  double p x n matrix (non-negative values)<br>
%         param: struct<br>
%            param.tol, param.max_it, param.verbose, param.L0,<br>
%            param.it0, param.lambda as for mexFistaFlat<br>
%            param.regul,  as for mexProximalFlat<br>
%<br>
% Output:<br>
%         W: double p x n matrix (non-negative values)</span><br>
<br>
<span class="c002">function</span> [W optim] = solverPoisson(Y,X,W0,param)<br>
param.ista=true;<br>
param.intercept=false;<br>
param.pos=true;<br>
param.linesearch_mode=2;<br>
param.loss=<span class="c003">'poisson'</span>;<br>
tabdelta=<span class="c002">logspace</span>(0,<span class="c002">log10</span>(param.delta),-<span class="c002">log10</span>(param.delta));<br>
<span class="c002">for</span> delta = tabdelta<br>
   param2=param;<br>
   param2.delta=delta;<br>
   [W optim]=mexFistaFlat(Y,X,W0,param2);<br>
   W0=W;<br>
<span class="c002">end</span></td></tr>
</table><p>
The following piece of code illustrates how to use this function.
</p><table class="lstframe c011"><tr><td class="mouselstlisting"><span class="c002">rand</span>(<span class="c003">'seed'</span>,0);<br>
<span class="c002">randn</span>(<span class="c003">'seed'</span>,0);<br>
<span class="c002">format</span> compact;<br>
<span class="c002">randn</span>(<span class="c003">'seed'</span>,0);<br>
param.numThreads=-1; <span class="c001">% all cores (-1 by default)</span><br>
param.verbose=true;   <span class="c001">% verbosity, false by default</span><br>
param.lambda=0.05; <span class="c001">% regularization parameter</span><br>
param.it0=10;      <span class="c001">% frequency for duality gap computations</span><br>
param.max_it=200; <span class="c001">% maximum number of iterations</span><br>
param.L0=0.01;<br>
param.tol=1e-6;<br>
param.delta=1e-5;<br>
<br>
X=<span class="c002">rand</span>(1000,2000);<br>
X=mexNormalize(X);<br>
Y=<span class="c002">rand</span>(1000,1);<br>
Y=mexNormalize(Y);<br>
W0=<span class="c002">zeros</span>(<span class="c002">size</span>(X,2),<span class="c002">size</span>(Y,2));<br>
<span class="c001">% Regression experiments <br>
% 100 regression problems with the same design matrix X.</span><br>
<span class="c002">fprintf</span>(<span class="c003">'\nVarious regression experiments\n'</span>);<br>
<span class="c002">fprintf</span>(<span class="c003">'\nFISTA + Regression l1\n'</span>);<br>
param.regul=<span class="c003">'l1'</span>;<br>
<span class="c002">tic</span><br>
[W optim_info]=solverPoisson(Y,X,W0,param);<br>
t=<span class="c002">toc</span>;<br>
<span class="c002">fprintf</span>(<span class="c003">'mean loss: %f, mean relative duality_gap: %f, time: %f, number of iterations: %f\n'</span>,<span class="c002">mean</span>(optim_info(1,:)),<span class="c002">mean</span>(optim_info(3,:)),t,<span class="c002">mean</span>(optim_info(4,:)));</td></tr>
</table>
<h3 class="subsection" id="sec42">5.13  Function mexIncrementalProx</h3>
<p>
This implements the incremental solver MISO [<a href="doc_spams010.html#mairal16">25</a>].
</p><table class="lstframe c011"><tr><td class="mouselstlisting"><span class="c001">% <br>
% Usage:  [W [optim]]=mexIncrementalProx(y,X,W0,param);<br>
%<br>
% Name: mexIncrementalProx<br>
%<br>
% Description: mexIncremrentalProx implements the incremental algorithm MISO<br>
% for composite optimization in a large scale setting.<br>
%        X is a design matrix of size p x n<br>
%        y is a vector of size n <br>
% WARNING, X is transposed compared to the functions mexFista*, and y is a vector<br>
%        param.lambda is a vector that contains nlambda different values of the<br>
%        regularization parameter (it can be a scalar, in that case nlambda=1)<br>
%        W0: is a dense matrix of size p x nlambda   It is in fact ineffective<br>
%        in the current release<br>
%        W: is the output, dense matrix of size p x nlambda<br>
%<br>
%         - if param.loss='square' and param.regul corresponds to a regularization<br>
%           function (currently 'l1' or 'l2'), the following problem is solved<br>
%           w = argmin (1/n)sum_{i=1}^n 0.5(y_i- x_i^T w)^2 + lambda psi(w)<br>
%         - if param.loss='logistic' and param.regul corresponds to a regularization<br>
%           function (currently 'l1' or 'l2'), the following problem is solved<br>
%           w = argmin (1/n)sum_{i=1}^n log(1+ exp(-y_ix_i^T w)) + lambda psi(w)<br>
%           Note that here, the y_i's should be -1 or +1 <br>
%          <br>
%         The current release does not handle intercepts<br>
%<br>
% Inputs: y: double dense vector of size n<br>
%         X: dense or sparse matrix of size p x n<br>
%         W0: dense matrix of size p x nlambda<br>
%         param: struct<br>
%           param.loss (choice of loss)<br>
%           param.regul (choice of regularization function)<br>
%           param.lambda : vector of size nlambda<br>
%           param.epochs: number of passes over the data<br>
%           param.minibatches: size of the mini-batches: recommended value is 1<br>
%              if X is dense, and min(n,ceil(1/density)) if X is sparse<br>
%           param.warm_restart : (path-following strategy, very efficient when<br>
%              providing an array of lambdas, false by default)<br>
%           param.normalized : (optional, can be set to true if the x_i's have<br>
%              unit l2-norm, false by default)<br>
%           param.strategy (optional, 3 by default)<br>
%                          0: no heuristics, slow  (only for comparison purposes)<br>
%                          1: adjust the constant L on 5% of the data <br>
%                          2: adjust the constant L on 5% of the data + unstable heuristics (this strategy does not work)<br>
%                          3: adjust the constant L on 5% of the data + stable heuristic (this is by far the best choice)<br>
%           param.numThreads (optional, number of threads)<br>
%           param.verbose (optional)<br>
%           param.seed (optional, choice of the random seed)<br>
%<br>
% Output:  W:  double dense p x nlambda matrix<br>
%          optim: optional, double dense 3 x nlambda matrix.<br>
%              first row: values of the objective functions.<br>
%              third row: computational time <br>
%<br>
% Author: Julien Mairal, 2013</span></td></tr>
</table><p>
The following piece of code illustrates how to use this function.
</p><table class="lstframe c011"><tr><td class="mouselstlisting"><span class="c002">clear</span> <span class="c002">all</span>;<br>
n=400000;<br>
p=1000;<br>
density=0.01;<br>
<br>
<span class="c001">% generate random data</span><br>
<span class="c002">format</span> compact;<br>
<span class="c002">randn</span>(<span class="c003">'seed'</span>,0);<br>
<span class="c002">rand</span>(<span class="c003">'seed'</span>,0);<br>
X=<span class="c002">sprandn</span>(p,n,density);<br>
mean_nrm=<span class="c002">mean</span>(<span class="c002">sqrt</span>(<span class="c002">sum</span>(X.^2)));<br>
X=X/mean_nrm;<br>
<br>
<span class="c001">% generate some true model</span><br>
z=double(<span class="c002">sign</span>(<span class="c002">full</span>(<span class="c002">sprandn</span>(p,1,0.05))));<br>
y=X<span class="c003">'*z;<br>
<br>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<br>
% EXPERIMENT 1: Lasso<br>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<br>
fprintf('</span>EXPERIMENT FOR LASSO\n<span class="c003">');<br>
nrm=sqrt(sum(y.^2));<br>
y=y+0.01*nrm*randn(n,1);    % add noise to the model<br>
nrm=sqrt(sum(y.^2));<br>
y=y*(sqrt(n)/nrm);<br>
<br>
% set optimization parameters<br>
clear param;<br>
param.regul='</span>l1<span class="c003">';        % many other regularization functions are available<br>
param.loss='</span>square<span class="c003">';     % only square and log are available<br>
param.numThreads=-1;    % uses all possible cores<br>
param.normalized=false;  % if the columns of X have unit norm, set to true.<br>
param.strategy=1;        % MISO with all heuristics<br>
                         % 0: no heuristics, slow  (only for comparison purposes)<br>
                         % 1: MISO1: adjust the constant L on 5% of the data  (good for non-strongly convex problems)<br>
                         % 2: adjust the constant L on 5% of the data + unstable heuristics (this strategy does not work)<br>
                         % 3: MISO2: adjust the constant L on 5% of the data + stable heuristic good for strongly-convex problems<br>
                         % 4: MISOmu: best for l2 regularization<br>
param.verbose=true;<br>
<br>
% set grid of lambda<br>
max_lambda=max(abs(X*y))/n;<br>
tablambda=max_lambda*(2^(1/4)).^(0:-1:-20);  % order from large to small<br>
tabepochs=[1 2 3 5 10];  % in this script, we compare the results obtained when changing the number of passes on the data.<br>
param.lambda=tablambda;<br>
%% The problem which will be solved is<br>
%%   min_beta  1/(2n) ||y-X'</span> <span class="c002">beta</span>||_2^2 + lambda ||<span class="c002">beta</span>||_1<br>
<span class="c001">% the problems for different lambdas are solve INDEPENDENTLY in parallel</span><br>
<span class="c002">fprintf</span>(<span class="c003">'EXPERIMENT: ALL LAMBDAS WITHOUT WARM RESTART\n'</span>);<br>
param.warm_restart=false;<br>
param.verbose=false;<br>
param.strategy=1;<br>
param.minibatches=<span class="c002">min</span>(n,<span class="c002">ceil</span>(1/density));  <span class="c001">% size of the minibatches, requires to store twice the size of X  with strategy=1</span><br>
obj=[];<br>
spar=[];<br>
<span class="c002">for</span> ii=1:<span class="c002">length</span>(tabepochs)<br>
   param.epochs=tabepochs(ii);   <span class="c001">% one pass over the data</span><br>
   <span class="c002">fprintf</span>(<span class="c003">'EXP WITH %d PASS OVER THE DATA\n'</span>,tabepochs(ii));<br>
   nlambdas=<span class="c002">length</span>(param.lambda);<br>
   Beta0=<span class="c002">zeros</span>(p,nlambdas);<br>
   <span class="c002">tic</span><br>
   [Beta tmp]=mexIncrementalProx(y,X,Beta0,param);<br>
   <span class="c002">toc</span><br>
   <span class="c002">fprintf</span>(<span class="c003">'Objective functions: \n'</span>);<br>
   obj=[obj; tmp(1,:)];<br>
   obj<br>
   <span class="c002">fprintf</span>(<span class="c003">'Sparsity: \n'</span>);<br>
   spar=[spar; <span class="c002">sum</span>(Beta ~= 0)];<br>
   spar<br>
<span class="c002">end<br>
<br>
<span class="c001">% the problems are here solved sequentially with warm restart<br>
% this seems to be the prefered choice.</span><br>
fprintf</span>(<span class="c003">'EXPERIMENT: SEQUENTIAL LAMBDAS WITH WARM RESTART\n'</span>);<br>
<span class="c002">fprintf</span>(<span class="c003">'A SINGLE CORE IS USED\n'</span>);<br>
param.warm_restart=true;<br>
param.num_threads=1;<br>
obj=[];<br>
spar=[];<br>
<span class="c002">for</span> ii=1:<span class="c002">length</span>(tabepochs)<br>
   param.epochs=tabepochs(ii);   <span class="c001">% one pass over the data</span><br>
   <span class="c002">fprintf</span>(<span class="c003">'EXP WITH %d PASS\n'</span>,tabepochs(ii));<br>
   nlambdas=<span class="c002">length</span>(param.lambda);<br>
   Beta0=<span class="c002">zeros</span>(p,nlambdas);<br>
   <span class="c002">tic</span><br>
   [Beta tmp]=mexIncrementalProx(y,X,Beta0,param);<br>
   <span class="c002">toc</span><br>
   <span class="c002">fprintf</span>(<span class="c003">'Objective functions: \n'</span>);<br>
   obj=[obj; tmp(1,:)];<br>
   obj<br>
   <span class="c002">fprintf</span>(<span class="c003">'Sparsity: \n'</span>);<br>
   spar=[spar; <span class="c002">sum</span>(Beta ~= 0)];<br>
   spar<br>
<span class="c002">end<br>
<br>
<span class="c001">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<br>
% EXPERIMENT 2: L2 logistic regression <br>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span><br>
fprintf</span>(<span class="c003">'EXPERIMENT FOR LOGISTIC REGRESSION + l2\n'</span>);<br>
y=<span class="c002">sign</span>(y);<br>
param.regul=<span class="c003">'l2'</span>;        <span class="c001">% many other regularization functions are available</span><br>
param.loss=<span class="c003">'logistic'</span>;     <span class="c001">% only square and log are available</span><br>
param.num_threads=-1;    <span class="c001">% uses all possible cores</span><br>
param.strategy=4;<br>
param.minibatches=1; <span class="c001">% with strategy=4, minibatches is better set to 1<br>
%param.strategy=3;        % MISO with all heuristics</span><br>
<span class="c002">fprintf</span>(<span class="c003">'EXPERIMENT: ALL LAMBDAS WITHOUT WARM RESTART\n'</span>);<br>
param.warm_restart=false;<br>
param.verbose=false;<br>
obj=[];<br>
spar=[];<br>
<span class="c002">for</span> ii=1:<span class="c002">length</span>(tabepochs)<br>
   param.epochs=tabepochs(ii);   <span class="c001">% one pass over the data</span><br>
   <span class="c002">fprintf</span>(<span class="c003">'EXP WITH %d PASS\n'</span>,tabepochs(ii));<br>
   nlambdas=<span class="c002">length</span>(param.lambda);<br>
   Beta0=<span class="c002">zeros</span>(p,nlambdas);<br>
   <span class="c002">tic</span><br>
   [Beta tmp]=mexIncrementalProx(y,X,Beta0,param);<br>
   <span class="c002">toc</span><br>
   yR=repmat(y,[1 nlambdas]);<br>
   <span class="c002">fprintf</span>(<span class="c003">'Objective functions: \n'</span>);<br>
   obj=[obj;tmp(1,:)];<br>
   obj<br>
<span class="c002">end<br>
<br>
<br>
fprintf</span>(<span class="c003">'EXPERIMENT FOR LOGISTIC REGRESSION + l2\n'</span>);<br>
y=<span class="c002">sign</span>(y);<br>
param.regul=<span class="c003">'l2'</span>;        <span class="c001">% many other regularization functions are available</span><br>
param.loss=<span class="c003">'logistic'</span>;     <span class="c001">% only square and log are available</span><br>
param.num_threads=-1;    <span class="c001">% uses all possible cores</span><br>
param.strategy=3;         <span class="c001">% MISO2</span><br>
param.minibatches=<span class="c002">min</span>(n,<span class="c002">ceil</span>(1/density));  <span class="c001">% size of the minibatches, requires to store twice the size of X </span><br>
<span class="c002">fprintf</span>(<span class="c003">'EXPERIMENT: ALL LAMBDAS WITHOUT WARM RESTART\n'</span>);<br>
param.warm_restart=false;<br>
param.verbose=false;<br>
obj=[];<br>
spar=[];<br>
<span class="c002">for</span> ii=1:<span class="c002">length</span>(tabepochs)<br>
   param.epochs=tabepochs(ii);   <span class="c001">% one pass over the data</span><br>
   <span class="c002">fprintf</span>(<span class="c003">'EXP WITH %d PASS\n'</span>,tabepochs(ii));<br>
   nlambdas=<span class="c002">length</span>(param.lambda);<br>
   Beta0=<span class="c002">zeros</span>(p,nlambdas);<br>
   <span class="c002">tic</span><br>
   [Beta tmp]=mexIncrementalProx(y,X,Beta0,param);<br>
   <span class="c002">toc</span><br>
   yR=repmat(y,[1 nlambdas]);<br>
   <span class="c002">fprintf</span>(<span class="c003">'Objective functions: \n'</span>);<br>
   obj=[obj;tmp(1,:)];<br>
   obj<br>
<span class="c002">end</span><br>
<br>
<span class="c001">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<br>
% EXPERIMENT 3: L1 logistic regression <br>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span><br>
<br>
y=<span class="c002">sign</span>(y);<br>
<span class="c002">fprintf</span>(<span class="c003">'EXPERIMENT FOR LOGISTIC REGRESSION + l1\n'</span>);<br>
param.regul=<span class="c003">'l1'</span>;        <span class="c001">% many other regularization functions are available</span><br>
param.loss=<span class="c003">'logistic'</span>;     <span class="c001">% only square and log are available</span><br>
param.num_threads=-1;    <span class="c001">% uses all possible cores</span><br>
param.strategy=1;<br>
param.minibatches=<span class="c002">min</span>(n,<span class="c002">ceil</span>(1/density));  <span class="c001">% size of the minibatches, requires to store twice the size of X </span><br>
<span class="c002">fprintf</span>(<span class="c003">'EXPERIMENT: ALL LAMBDAS WITHOUT WARM RESTART\n'</span>);<br>
param.warm_restart=false;<br>
param.verbose=false;<br>
param.lambda=tablambda;<br>
obj=[];<br>
spar=[];<br>
<span class="c002">for</span> ii=1:<span class="c002">length</span>(tabepochs)<br>
   param.epochs=tabepochs(ii);   <span class="c001">% one pass over the data</span><br>
   <span class="c002">fprintf</span>(<span class="c003">'EXP WITH %d PASS\n'</span>,tabepochs(ii));<br>
   nlambdas=<span class="c002">length</span>(param.lambda);<br>
   Beta0=<span class="c002">zeros</span>(p,nlambdas);<br>
   <span class="c002">tic</span><br>
   [Beta tmp]=mexIncrementalProx(y,X,Beta0,param);<br>
   <span class="c002">toc</span><br>
   yR=repmat(y,[1 nlambdas]);<br>
   <span class="c002">fprintf</span>(<span class="c003">'Objective functions: \n'</span>);<br>
   obj=[obj;tmp(1,:)];<br>
   obj<br>
<span class="c002">end</span></td></tr>
</table>
<h3 class="subsection" id="sec43">5.14  Function mexStochasticProx</h3>
<p>
This implements the stochastic proximal gradient solver [<a href="doc_spams010.html#mairal15">26</a>].
</p><table class="lstframe c011"><tr><td class="mouselstlisting"><span class="c001">% <br>
% Usage:  [W [W2]]=mexStochasticProx(y,X,W0,param);<br>
%<br>
% Name: mexStochasticProx<br>
%<br>
% Description: mexStochasticProx implements a proximal MM stochastic algorithm <br>
% for composite optimization in a large scale setting.<br>
%        X is a design matrix of size p x n<br>
%        y is a vector of size n <br>
% WARNING, X is transposed compared to the functions mexFista*, and y is a vector<br>
%        param.lambda is a vector that contains nlambda different values of the<br>
%        regularization parameter (it can be a scalar, in that case nlambda=1)<br>
%        W0: is a dense matrix of size p x nlambda   It is in fact ineffective<br>
%        in the current release<br>
%        W: is the output, dense matrix of size p x nlambda<br>
%<br>
%         - if param.loss='square' and param.regul corresponds to a regularization<br>
%           function (currently 'l1' or 'l2'), the following problem is solved<br>
%           w = argmin (1/n)sum_{i=1}^n 0.5(y_i- x_i^T w)^2 + lambda psi(w)<br>
%         - if param.loss='logistic' and param.regul corresponds to a regularization<br>
%           function (currently 'l1' or 'l2'), the following problem is solved<br>
%           w = argmin (1/n)sum_{i=1}^n log(1+ exp(-y_ix_i^T w)) + lambda psi(w)<br>
%           Note that here, the y_i's should be -1 or +1 <br>
%          <br>
%         The current release does not handle intercepts<br>
%<br>
% Inputs: y: double dense vector of size n<br>
%         X: dense or sparse matrix of size p x n<br>
%         W0: dense matrix of size p x nlambda<br>
%         param: struct<br>
%           param.loss (choice of loss)<br>
%           param.regul (choice of regularization function)<br>
%           param.lambda : vector of size nlambda<br>
%           param.iters : number of iterations (n corresponds to one pass over the data)<br>
%           param.minibatches: size of the mini-batches: recommended value is 1<br>
%           param.normalized : (optional, can be set to true if the x_i's have<br>
%              unit l2-norm, false by default)<br>
%           param.weighting_mode : (optional, 1 by default),<br>
%                0:  w_t = (t_0+1)/(t+t_0)<br>
%                1:  w_t = ((t_0+1)/(t+t_0))^(0.75)<br>
%                2:  w_t = ((t_0+1)/(t+t_0))^(5)<br>
%           param.averaging_mode: (optional, false by default)<br>
%                0: no averaging<br>
%                1: first averaging mode for W2<br>
%                2: second averaging mode for W2<br>
%                WARNING: averaging can be very slow for sparse solutions<br>
%           param.determineEta (optional, automatically choose the parameters of the<br>
%             learning weights w_t, true by default) <br>
%           param.t0 (optional, set up t0 for weights w_t = ((1+t0)/(t+t0))^(alpha)<br>
%           param.numThreads (optional, number of threads)<br>
%           param.verbose (optional)<br>
%           param.seed (optional, choice of the random seed)<br>
%<br>
% Output:  W:  double dense p x nlambda matrix (contains the solution without averaging)<br>
%          W2:  double dense p x nlambda matrix (contains the solution with averaging)<br>
%<br>
% Author: Julien Mairal, 2013</span></td></tr>
</table><p>
The following piece of code illustrates how to use this function.
</p><table class="lstframe c011"><tr><td class="mouselstlisting"><span class="c002">clear</span> <span class="c002">all</span>;<br>
n=400000;<br>
p=1000;<br>
density=0.01;<br>
<br>
<span class="c001">% generate random data</span><br>
<span class="c002">format</span> compact;<br>
<span class="c002">randn</span>(<span class="c003">'seed'</span>,0);<br>
<span class="c002">rand</span>(<span class="c003">'seed'</span>,0);<br>
X=<span class="c002">sprandn</span>(p,n,density);<br>
nrm=<span class="c002">sqrt</span>(<span class="c002">sum</span>(X.^2));<br>
X=X/<span class="c002">mean</span>(nrm);<br>
<span class="c001">%mean_nrm=mean();<br>
%X=X/mean_nrm;<br>
<br>
% generate some true model</span><br>
z=double(<span class="c002">sign</span>(<span class="c002">full</span>(<span class="c002">sprandn</span>(p,1,0.05))));<br>
y=X<span class="c003">'*z;<br>
<br>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<br>
% EXPERIMENT 1: Lasso<br>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<br>
fprintf('</span>EXPERIMENT FOR LASSO\n<span class="c003">');<br>
nrm=sqrt(sum(y.^2));<br>
y=y+0.01*nrm*randn(n,1);    % add noise to the model<br>
nrm=sqrt(sum(y.^2));<br>
y=y*(sqrt(n)/nrm);<br>
<br>
clear param;<br>
param.regul='</span>l1<span class="c003">';        % many other regularization functions are available<br>
param.loss='</span>square<span class="c003">';     % only square and log are available<br>
param.numThreads=-1;    % uses all possible cores<br>
param.normalized=false;  % if the columns of X have unit norm, set to true.<br>
param.averaging_mode=0;  % no averaging, averaging was not really useful in experiments<br>
param.weighting_mode=0;  % weights are in O(1/sqrt(n))  (see help mexStochasticProx)<br>
param.optimized_solver=true;  % best is not to touch this option<br>
param.verbose=false;<br>
<br>
% set grid of lambda<br>
max_lambda=max(abs(X*y))/n;<br>
tablambda=max_lambda*(2^(1/4)).^(0:-1:-20);  % order from large to small<br>
param.lambda=tablambda;    % best to order from large to small<br>
tabepochs=[1 2 3 5 10];  % in this script, we compare the results obtained when varying the number of passes over the data.<br>
<br>
%% The problem which will be solved is<br>
%%   min_beta  1/(2n) ||y-X'</span> <span class="c002">beta</span>||_2^2 + lambda ||<span class="c002">beta</span>||_1<br>
<span class="c002">fprintf</span>(<span class="c003">'EXPERIMENT: ALL LAMBDAS IN PARALLEL\n'</span>);<br>
<span class="c001">% we try different experiments when varying the number of epochs.<br>
% the problems for different lambdas are solve INDEPENDENTLY in parallel</span><br>
obj=[];<br>
objav=[];<br>
<span class="c002">for</span> ii=1:<span class="c002">length</span>(tabepochs)<br>
   param.iters=tabepochs(ii)*n;<br>
   <span class="c002">fprintf</span>(<span class="c003">'EXP WITH %d PASS\n'</span>,tabepochs(ii));<br>
   nlambdas=<span class="c002">length</span>(param.lambda);<br>
   Beta0=<span class="c002">zeros</span>(p,nlambdas);<br>
   <span class="c002">tic</span><br>
   [Beta tmp]=mexStochasticProx(y,X,Beta0,param);<br>
   <span class="c002">toc</span><br>
   yR=repmat(y,[1 nlambdas]);<br>
   <span class="c002">fprintf</span>(<span class="c003">'Objective functions: \n'</span>);<br>
   obj=[obj; 0.5*<span class="c002">sum</span>((yR-X<span class="c003">'*Beta).^2)/n+param.lambda.*sum(abs(Beta))];<br>
   obj<br>
   if param.averaging_mode<br>
      objav=[objav; 0.5*sum((yR-X'</span>*tmp).^2)/n+param.lambda.*<span class="c002">sum</span>(<span class="c002">abs</span>(tmp))];<br>
      objav<br>
   <span class="c002">end</span><br>
   <span class="c002">fprintf</span>(<span class="c003">'Sparsity: \n'</span>);<br>
   <span class="c002">sum</span>(Beta ~= 0)<br>
<span class="c002">end<br>
<br>
<span class="c001">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<br>
% EXPERIMENT 2: L2 logistic regression <br>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span><br>
fprintf</span>(<span class="c003">'EXPERIMENT FOR LOGISTIC REGRESSION + l2\n'</span>);<br>
y=<span class="c002">sign</span>(y);<br>
param.regul=<span class="c003">'l2'</span>;<br>
param.loss=<span class="c003">'logistic'</span>;<br>
obj=[];<br>
objav=[];<br>
<span class="c002">for</span> ii=1:<span class="c002">length</span>(tabepochs)<br>
   param.iters=tabepochs(ii)*n;   <span class="c001">% one pass over the data</span><br>
   <span class="c002">fprintf</span>(<span class="c003">'EXP WITH %d PASS\n'</span>,tabepochs(ii));<br>
   nlambdas=<span class="c002">length</span>(param.lambda);<br>
   Beta0=<span class="c002">zeros</span>(p,nlambdas);<br>
   <span class="c002">tic</span><br>
   [Beta tmp]=mexStochasticProx(y,X,Beta0,param);<br>
   <span class="c002">toc</span><br>
   yR=repmat(y,[1 nlambdas]);<br>
   <span class="c002">fprintf</span>(<span class="c003">'Objective functions: \n'</span>);<br>
   obj=[obj;<span class="c002">sum</span>(<span class="c002">log</span>(1.0+<span class="c002">exp</span>(-yR .* (X<span class="c003">'*Beta))))/n+0.5*param.lambda.*sum(abs(Beta.^2))];<br>
   obj<br>
   if param.averaging_mode<br>
      objav=[objav; sum(log(1.0+exp(-yR .* (X'</span>*tmp))))/n+0.5*param.lambda.*<span class="c002">sum</span>(<span class="c002">abs</span>(tmp))];<br>
      objav<br>
   <span class="c002">end<br>
end<br>
<br>
<span class="c001">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<br>
% EXPERIMENT 3: L1 logistic regression <br>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span><br>
fprintf</span>(<span class="c003">'EXPERIMENT FOR LOGISTIC REGRESSION + l1\n'</span>);<br>
y=<span class="c002">sign</span>(y);<br>
param.regul=<span class="c003">'l1'</span>;        <span class="c001">% many other regularization functions are available</span><br>
param.loss=<span class="c003">'logistic'</span>;     <span class="c001">% only square and log are available</span><br>
obj=[];<br>
objav=[];<br>
<span class="c002">for</span> ii=1:<span class="c002">length</span>(tabepochs)<br>
   param.iters=tabepochs(ii)*n;   <span class="c001">% one pass over the data</span><br>
   <span class="c002">fprintf</span>(<span class="c003">'EXP WITH %d PASS\n'</span>,tabepochs(ii));<br>
   nlambdas=<span class="c002">length</span>(param.lambda);<br>
   Beta0=<span class="c002">zeros</span>(p,nlambdas);<br>
   <span class="c002">tic</span><br>
   [Beta tmp]=mexStochasticProx(y,X,Beta0,param);<br>
   <span class="c002">toc</span><br>
   yR=repmat(y,[1 nlambdas]);<br>
   <span class="c002">fprintf</span>(<span class="c003">'Objective functions: \n'</span>);<br>
   obj=[obj; <span class="c002">sum</span>(<span class="c002">log</span>(1.0+<span class="c002">exp</span>(-yR .* (X<span class="c003">'*Beta))))/n+param.lambda.*sum(abs(Beta))];<br>
   obj<br>
   if param.averaging_mode<br>
      objav=[objav; sum(log(1.0+exp(-yR .* (X'</span>*tmp))))/n+param.lambda.*<span class="c002">sum</span>(<span class="c002">abs</span>(tmp))];<br>
      objav<br>
   <span class="c002">end</span><br>
   <span class="c002">fprintf</span>(<span class="c003">'Sparsity: \n'</span>);<br>
   <span class="c002">sum</span>(Beta ~= 0)<br>
<span class="c002">end</span></td></tr>
</table>
<hr>
<a href="doc_spams005.html"><img src="previous_motif.gif" alt="Previous"></a>
<a href="index.html"><img src="contents_motif.gif" alt="Up"></a>
<a href="doc_spams007.html"><img src="next_motif.gif" alt="Next"></a>
</body>
</html>
